2014-12-01 3 views
2

У меня длинная последовательность последовательности ДНК, и мне нужно найти области, состоящие из двух палиндромных последовательностей, фланкирующих спейсерную последовательность.Поиск инвертированных повторов в последовательности ДНК

вход является:

cgtacacgagtagtcgtagctgtcagtcgatcgtacgtacgtagctgctgtagcactatcgaccccacacgtgtgtacacgatgcacagtcgtctatcacatgctagcgctgcccgtacgGATGGCCAAGGCCATCcgatcgctagctagcgccgcgcgtagcccgatcgagacatgctagcagttgtgctgatgtcgagatagctgtgatgcgatgctagcgccgcctagccgcctcgtgtaggctggatgcga tcgatcgatgctagcggcgcgatcga tgcactagccgtagcgctagctgatcgatcgtaGATGGCCAAGGCCATCcgcgtagatacgacatgccgggggtatataa

Это мой код:

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

с параметрами: $ мин = 6, $ макс = 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. Может ли кто-нибудь просветить меня, что случилось с моим кодом? Благодарю.

+0

Я попытался использовать другие последовательности ДНК, и код, казалось, работал, но кроме этого. – zebra

+0

Извините, я не могу понять, что делает ваш код, и что именно не работает. Не могли бы вы немного разъяснить? – Sobrique

+0

Извините, я не комментировал свой код. Я намерен идентифицировать регион с двумя фланкирующими палиндромами ДНК. Например, GAATTC является палиндром его комплементарной нити при чтении в обратном направлении. Пример региона может быть GAATTC ----- GAATTC или GGGC ---- GCCC, где - нуклеотид в спейсерной области. – zebra

ответ

1

Ваша проблема заключается в том, что ваше регулярное выражение ищет палиндромы в любом месте последовательности, а не только в месте смещения. «ATCGATCG» встречается со смещением 26, поэтому он соответствует. Вам нужно добавить некоторую позиционную информацию в ваше регулярное выражение. Попробуйте что-то вроде

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

Хотя обратите внимание, что это не будет работать в течение длинных последовательностей ДНК, где '$ offset' > 32766. Лучше всего было бы искать только в окне: '$ search_window = substr ($ dna, $ pos, 500);' if if $ $$$$$$$$$$$$$$$$$$$$$$$$$$ ACGT] {$ min_spacer, $ max_spacer}) ($ inverted_region)) /) {' – fugu

1

Вот одно решение, он использует экспериментальную (??{}) функцию, которая, как говорят, чтобы изменить в течение длительного времени, но не имеет еще комментариев.

Как это работает: он вызывает подпрограмму 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 

Кроме того, я не уверен, что это проблема, но вы могли бы производить больше palidrome последовательности, которые просто сдвигаются в распорки с этим решением:

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