如何使用 Perl 搜索字符串中生物信息学主题的不同变体？

问题描述 投票：0回答：3

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

perl string bioinformatics
3个回答
4
投票

我不知道我使用的

'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 


    



      

      
      
      

      

        

          

            

              
1
投票

              
            

          

        

        



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




      

      
      
      

      

        

          

            

              
0
投票

              
            

          

        

        



以下是如何使用它的示例：



#!/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. 

  2. 共识序列
    3. 

  3. 找到的单位数量
    4. 

  4. 质量指标（满分 1）
    5. 

  5. 重复单元之间用空格分隔
    6. 

  6. 


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



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



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

    



      

      
    

  

  

    
最新问题


  





  

    © www.soinside.com 2019 - 2024. All rights reserved.