在 DNA 序列中寻找反向重复答案

【问题标题】：Finding inverted repeats in DNA sequence在 DNA 序列中寻找反向重复
【发布时间】：2014-12-01 04:31:53
【问题描述】：

我有一长串 DNA 序列，我需要找到由两个回文序列组成的区域，它们位于间隔序列的两侧。

输入是：

cgtacacgagtagtcgtagctgtcagtcgatcgtacgtacgtagctgctgtagcactatcgaccccacacgtgtgtacacgatgcacagtcgtctatcacatgctagcgctgcccgtacgGATGGCCAAGGCCATCcgatcgctagctagcgccgcgcgtagcccgatcgagacatgctagcagttgtgctgatgtcgagatagctgtgatgcgatgctagcgccgcctagccgcctcgtgtaggctggatgcga的 tcgatcgatgctagcggcgcgatcga 强> tgcactagccgtagcgctagctgatcgatcgtaGATGGCCAAGGCCATCcgcgtagatacgacatgccgggggtatataa P>

这是我的代码：

use strict;
use warnings;

my $input= $ARGV[0];
chomp $input;
open (my $fh_in, "<", $input) or die "Cannot open file $input $!";
my $dna= <$fh_in>;
chomp $dna;

#######################################################################################

if ($dna=~ /[^ACGT]/gi) {
        print "This is not a valid DNA sequence, it has unknown base(s)\n";
}

$dna=~ tr/[acgt]/[ACGT]/;


######################################################################################

print "Minimum length of palindromic sequence?\n";
my $min= <STDIN>;
chomp $min;

print "Maximum length of palindromic sequence?\n";
my $max= <STDIN>;
chomp $max;

print "Minimum length of spacer region?\n";
my $min_spacer= <STDIN>;
chomp $min_spacer;

print "Maximum length of spacer region?\n";
my $max_spacer= <STDIN>;
chomp $max_spacer;

###################################################################################### 

my $dna_length= length($dna);
my ($length , $offset , $string_1 , $string_2);

for ($offset= 0 ; $offset <= $dna_length-$max-$max-$max_spacer ; $offset++) {
        for ($length= $min ; $length <= $max ; $length++) {
                $string_1= substr ($dna, $offset, $length);
                $string_2= reverse $string_1;
                $string_2=~ tr/[ACGT]/[TGCA]/;
                if ($dna=~ /(($string_1)([ACGT]{$min_spacer,$max_spacer})($string_2))/) {
                        print "IR starts at $offset => $2***$3***$4\n$1\n\n";
                }
        }
}

带参数： $min = 6 , $max = 12 , $min_spacer = 4 , $max_spacer = 12 我得到的输出是：

IR starts at 26 => TCGATCG***ATGCTAGCGGCG***CGATCGA
TCGATCGATGCTAGCGGCGCGATCGA

IR starts at 27 => CGATCG***ATGCTAGCGGCG***CGATCG
CGATCGATGCTAGCGGCGCGATCG

IR starts at 118 => CGGATG***GCCAAGGC***CATCCG
CGGATGGCCAAGGCCATCCG

IR starts at 118 => CGGATGG***CCAAGG***CCATCCG
CGGATGGCCAAGGCCATCCG

IR starts at 118 => CGGATGGC***CAAG***GCCATCCG
CGGATGGCCAAGGCCATCCG

IR starts at 119 => GGATGG***CCAAGG***CCATCC
GGATGGCCAAGGCCATCC

IR starts at 119 => GGATGGC***CAAG***GCCATCC
GGATGGCCAAGGCCATCC

IR starts at 120 => GATGGC***CAAG***GCCATC
GATGGCCAAGGCCATC

IR starts at 136 => CGATCG***ATGCTAGCGGCG***CGATCG
CGATCGATGCTAGCGGCGCGATCG

IR starts at 164 => CGATCG***ATGCTAGCGGCG***CGATCG
CGATCGATGCTAGCGGCGCGATCG

IR starts at 252 => CGATCG***ATGCTAGCGGCG***CGATCG
CGATCGATGCTAGCGGCGCGATCG

IR starts at 254 => ATCGAT***GCTAGCGGCGCG***ATCGAT
ATCGATGCTAGCGGCGCGATCGAT

IR starts at 254 => ATCGATCG***ATGCTAGCGGCG***CGATCGAT
ATCGATCGATGCTAGCGGCGCGATCGAT

IR starts at 255 => TCGATCG***ATGCTAGCGGCG***CGATCGA
TCGATCGATGCTAGCGGCGCGATCGA

IR starts at 256 => CGATCG***ATGCTAGCGGCG***CGATCG
CGATCGATGCTAGCGGCGCGATCG

IR starts at 258 => ATCGAT***GCTAGCGGCGCG***ATCGAT
ATCGATGCTAGCGGCGCGATCGAT

IR starts at 274 => CGATCG***ATGCTAGCGGCG***CGATCG
CGATCGATGCTAGCGGCGCGATCG

IR starts at 276 => ATCGAT***GCTAGCGGCGCG***ATCGAT
ATCGATGCTAGCGGCGCGATCGAT

IR starts at 304 => ATCGAT***GCTAGCGGCGCG***ATCGAT
ATCGATGCTAGCGGCGCGATCGAT

IR starts at 304 => ATCGATCG***ATGCTAGCGGCG***CGATCGAT
ATCGATCGATGCTAGCGGCGCGATCGAT

IR starts at 305 => TCGATCG***ATGCTAGCGGCG***CGATCGA
TCGATCGATGCTAGCGGCGCGATCGA

IR starts at 306 => CGATCG***ATGCTAGCGGCG***CGATCG
CGATCGATGCTAGCGGCGCGATCG

IR starts at 314 => GATGGC***CAAG***GCCATC
GATGGCCAAGGCCATC

但是，当我检查我的第一个命中的区域（在输入中以粗体突出显示）时，这个命中的偏移量似乎不在位置 26。谁能告诉我我的代码有什么问题？谢谢。

【问题讨论】：

我尝试使用其他 DNA 序列，代码似乎有效，但除了这个。
抱歉，我无法弄清楚您的代码在做什么，以及究竟是什么不起作用。你能澄清一下吗？
对不起，我没有评论我的代码。我打算确定一个具有两个侧翼 DNA 回文的区域。例如，当反向阅读时，GAATTC 是其互补链的回文。区域示例可以是 GAATTC-----GAATTC 或 GGGC----GCCC，其中 - 是间隔区中的核苷酸。
您没有得到的预期输出是什么？恐怕我 - 像该网站的大多数读者一样 - 发现代码很容易，而 DNA 很难 :)
你们为什么要使用 perl 正则表达式引擎来完成您的任务？用你大大减少的字母表对我来说这似乎是非常低效的。

标签： regex perl bioinformatics dna-sequence

【解决方案1】：

您的问题是您的正则表达式正在序列中的任何位置寻找回文，而不仅仅是在偏移的位置。 “ATCGATCG”出现在偏移量 26 处，因此匹配。您需要在正则表达式中添加一些位置信息。尝试类似

/^[ACGT]{$offset}(($string_1)([ACGT]{$min_spacer,$max_spacer})($string_2))/

【讨论】：

虽然请注意，这不适用于长 DNA 序列，其中 $offset 大于 32766。更好的方法可能是仅在窗口内搜索：$search_window = substr($dna, $pos, 500);if ($search_window =~ /^[ACGT](($region)([ACGT]{$min_spacer,$max_spacer})($inverted_region))/) {

【解决方案2】：

这里有一个解决方案，它使用了实验性的(??{}) 功能，据说已经改变了很长时间，但还没有。

它是如何工作的：它从正则表达式中调用子例程convert，并将第一个匹配组转换为输出字符串的所需正则表达式。其余的（回溯等）由正则表达式引擎处理。遗憾的是，将变量作为定界长度进行插值并不适合正则表达式解析，因此我不得不使用字符串来做到这一点。如果可能，请不要这样做。

use warnings;
use strict;
use 5.01;
use re 'eval'; # needed, because of (??{})

my %c=
  (min_pali   =>    (shift) // 6,
   max_pali   =>    (shift) // 12,
   min_spacer =>    (shift) // 4,
   max_spacer =>    (shift) // 12,
  );

my $re1 = "(.{$c{min_pali},$c{max_pali}})(.{$c{min_spacer},$c{max_spacer}})(??{convert})";
while(<DATA>){
  chomp;
  $_ = uc $_;
  my $converted;
  sub convert {
    my $var = reverse $1;
    $var =~ tr{ACGT}{TGCA};
    $converted =  $var;
  }
  while (/$re1/g) {
    printf "%3d => %s**%s**%s\n", $-[0],$1,$2,$converted;
    pos = $-[0] + 1; # start next match one character after the last match start
  }
}

__DATA__
cgtacacgagtagtcgtagctgtcagtcgatcgtacgtacgtagctgctgtagcactatcgaccccacacgtgtgtacacgatgcacagtcgtctatcacatgctagcgctgcccgtacgGATGGCCAAGGCCATCcgatcgctagctagcgccgcgcgtagcccgatcgagacatgctagcagttgtgctgatgtcgagatagctgtgatgcgatgctagcgccgcctagccgcctcgtgtaggctggatgcgatcgatcgatgctagcggcgcgatcgatgcactagccgtagcgctagctgatcgatcgtaGATGGCCAAGGCCATCcgcgtagatacgacatgccgggggtatataa

输出：

118 => CGGATG**GCCAAGGC**CATCCG
119 => GGATGG**CCAAGG**CCATCC
120 => GATGGC**CAAG**GCCATC
254 => ATCGATCG**ATGCTAGCGGCG**CGATCGAT
255 => TCGATCG**ATGCTAGCGGCG**CGATCGA
256 => CGATCG**ATGCTAGCGGCG**CGATCG
258 => ATCGAT**GCTAGCGGCGCG**ATCGAT
314 => GATGGC**CAAG**GCCATC

我也不确定这是否是个问题，但您可以使用此解决方案生成更长的回文序列，然后将其转移到间隔中：

 Assuming length 2 – 4, spacer= 2 – 4 (X's are unintresting bits)
 ACACAXXTGTGT => ACAC**AXXT**GTGT

【讨论】：