2

我有一个带有不同变体的串联重复的程序输出。是否可以（在字符串中）搜索主题并告诉程序找到最大“3”不匹配/插入/删除的所有变体？

4

3 回答 3

4

我将利用提供的非常有限的信息来解决这个问题。

首先，简短友好的社论：

<editorial>

请学习如何提出一个好问题以及如何准确。

请至少：

  • 在不提供链接或精确定义的情况下避免使用特定领域的术语，例如“基序”、“串联重复”和“碱基对”；
  • 说出目标是什么，以及到目前为止你做了什么；
  • 重要提示：提供输入和所需输出的清晰示例。

这对 SO 的潜在帮助者没有帮助，必须在评论中回答 20 个问题才能尝试理解您的问题！我花了更多时间试图弄清楚你在问什么，而不是回答它。

</editorial>

以下程序在 1,000 个元素长的数组中生成 2 个字符对 5,428 对长的字符串。我意识到您更有可能从文件中读取这些内容，但这只是一个示例。显然，您将用来自任何来源的实际数据替换随机字符串。

我不知道'AT','CG','TC','CA','TG','GC','GG'我使用的是否是合法的碱基对组合。（我睡过生物学...）如果您想生成合法的随机字符串进行测试，只需将map块对编辑为合法对并将 更改为对数。7

如果该offset点的子字符串是 3 个或更少的差异，则数组元素（标量值）存储在哈希值部分的匿名数组中。哈希的关键部分是近似匹配的子字符串。除了数组元素，这些值可以是文件名、Perl 数据引用或您希望与您的主题相关联的其他相关引用。

虽然我刚刚查看了字符串之间逐个字符的差异，但您可以通过将行替换为foreach my $j (0..$#a1) { $diffs++ unless ($a1[$j] eq $a2[$j]); }适用于您的问题的比较逻辑来放置您需要查看的任何特定逻辑。我不知道mismatches/insertions/deletions在你的字符串中是如何表示的，所以我把它作为练习留给读者。也许来自 CPAN 的Algorithm::Diff或String::Diff ？

很容易修改这个程序，让键盘输入$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
于 2010-09-18T17:38:43.183 回答
1

我相信在 BioPerl 中有这类事情的例程。

无论如何，如果您在生物信息学堆栈交换 BioStar 上问过这个问题，您可能会得到更好的答案。

于 2010-09-19T15:11:21.653 回答
0

当我刚开始学习 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

如您所见，它允许插入缺失和 SNP。这些列按顺序排列：

  1. 起始位置
  2. 停止位置
  3. 共识序列
  4. 找到的单位数
  5. 质量指标 1
  6. 由空格分隔的重复单元

请注意，很容易提供将输出垃圾/无关紧要的“重复”的参数（正如您从上面的输出中看到的那样），但是如果您知道如何提供好的参数，它可以找到您在找到时设置的内容。

不幸的是，该软件包不公开。我从来没有费心让它可用，因为它太慢了，甚至不适合原核大小的基因组搜索（尽管它对单个基因是可行的）。在我刚开始编码的日子里，我开始添加一个功能来将“状态”作为输入，这样我就可以在序列的各个部分并行运行它，但我从来没有完成过，一旦我学会了哈希会让它变得更快。到那时，我已经转向其他项目。但是，如果它适合您的需要，请给我留言，我可以通过电子邮件将副本发送给您。

它只有不到 1000 行代码，但它有很多花里胡哨，例如允许 IUPAC 歧义代码 (BDHVRYKMSWN)。它适用于氨基酸和核酸。它过滤掉内部重复（例如，不将 TTTT 或 ATAT 报告为 4nt 共识）。

于 2016-03-31T21:56:10.907 回答