9 de julio de 2014

regenerando DNA degenerado

Hola,
durante la reciente visita de mi colega Pablo Vinuesa al laboratorio pasamos ratos escribiendo código en Perl con el fin de diseñar de manera automática, para grandes conjuntos de secuencias de DNA previamente alineadas, parejas de primers (cebadores de PCR) adecuadas para estudios de microbiología ambiental, siguiendo principios bien conocidos ya. En cuanto probamos el código con unos cuantos ejemplos observamos que a menudo los primers para algunas regiones de los alineamientos múltiples eran degenerados. Como se muestra en la figura, eso significa que algunas posiciones de los cebadores no podían ser nucleótidos únicos si no que tenían que ser combinaciones de 2 o más para poder hibridarse con las secuencias utlizadas para su diseño.


Pareja de primers degenerados que definen un amplicón a partir de secuencias de DNA alineadas. De acuerdo con la nomenclatura IUPAC, el de la izquierda (fwd) puede escribirse como GAYTST y el de la derecha (rev) como GHGKAG. Figura tomada de http://acgt.cs.tau.ac.il/hyden.
A la hora de evaluar parejas de primers nos encontramos con que los degenerados en realidad tendrían que ser examinados por medio de cada una de las secuencias implícitas en el código degenerado. Y por tanto, tuvimos que buscar una manera de regenerar las secuencias de nucleótidos correspondientes a un primer degenerado, que es de lo que va esta entrada. El siguiente código en Perl hace precisamente esto:

 #!/usr/bin/perl  
 use strict;  
 my $degenerate_sequence = $ARGV[0] || die "# usage: $0 <degenerate DNA sequence>\n";  
 my $regenerated_seqs = regenerate($degenerate_sequence);  
 foreach my $seq (@$regenerated_seqs)  
 {  
    print "$seq\n";  
 }  
 # returns a ref to list of all DNA sequences implied in a degenerate oligo  
 # adapted from Amplicon (Simon Neil Jarman, http://sourceforge.net/projects/amplicon)  
 sub regenerate  
 {  
    my ($primerseq) = @_;  
    my %IUPACdegen = (   
    'A'=>['A'],'C'=>['C'], 'G'=>['G'], 'T'=>['T'],   
    'R'=>['A','G'], 'Y'=>['C','T'], 'S'=>['C','G'], 'W'=>['A','T'], 'K'=>['G','T'], 'M'=>['A','C'],  
    'B'=>['C','G','T'], 'D'=>['A','G','T'], 'V'=>['A','C','G'], 'H'=>['A','C','T'],   
    'N'=>['A','C','G','T']   
    );  
    my @oligo = split(//,uc($primerseq));   
    my @grow = ('');  
    my @newgrow = ('');  
    my ($res,$degen,$x,$n,$seq);  
    foreach $res (0 .. $#oligo){  
       $degen = $IUPACdegen{$oligo[$res]};   
       if($#{$degen}>0){  
          $x = 0;  
          @newgrow = @grow;  
          while($x<$#{$degen}){  
             push(@newgrow,@grow);  
             $x++;     
          }     
          @grow = @newgrow;  
       }  
       $n=$x=0;   
       foreach $seq (0 .. $#grow){  
          $grow[$seq] .= $degen->[$x];   
          $n++;  
          if($n == (scalar(@grow)/scalar(@$degen))){  
             $n=0;  
             $x++;  
          }     
       }  
    }  
    return \@grow;     
 }            

Podemos probarlo con las secuencias de la figura:

$ perl degenprimer.pl GAYTST
GACTCT
GATTCT
GACTGT
GATTGT

 y

$ perl degenprimer.pl GHGKAG
GAGGAG
GCGGAG
GTGGAG
GAGTAG
GCGTAG
GTGTAG

Hasta otra,
Bruno

1 comentario:

  1. Os dejo el mismo código en python ;)

    def disambiguate(seq,stype='dna'):
    '''Retrieves all possible sequences from a IUPAC sequence'''

    nt_iupac = dict()

    nt_iupac = {
    'A' : ['A'], 'a' : ['a'],
    'C' : ['C'], 'c' : ['c'],
    'G' : ['G'], 'g' : ['g'],
    'T' : ['T'], 't' : ['t'],
    'U' : ['U'], 'u' : ['u'],
    'M' : ['A', 'C'], 'm' : ['a', 'c'],
    'R' : ['A', 'G'], 'r' : ['a', 'g'],
    'W' : ['A', 'T'], 'w' : ['a', 't'],
    'S' : ['C', 'G'], 's' : ['c', 'g'],
    'Y' : ['C', 'T'], 'y' : ['c', 't'],
    'K' : ['G', 'T'], 'k' : ['g', 't'],
    'V' : ['A', 'C', 'G'], 'v' : ['a', 'c', 'g'],
    'H' : ['A', 'C', 'T'], 'h' : ['a', 'c', 't'],
    'D' : ['A', 'G', 'T'], 'd' : ['a', 'g', 't'],
    'B' : ['C', 'G', 'T'], 'b' : ['c', 'g', 't'],
    'N' : ['A', 'C', 'G', 'T'], 'n' : ['a', 'c', 'g', 't'],
    'X' : ['A', 'C', 'G', 'T'], 'x' : ['a', 'c', 'g', 't']
    }

    count_combs = number_combs(seq)
    index_combs = count_combs
    comb_seqs = ['']*count_combs
    for nt in seq:
    comb_nts = nt_iupac[nt]
    index_combs = index_combs/len(comb_nts)
    count_nt_combs = 0
    index = 0
    for i in range(count_combs):
    count_nt_combs += 1
    comb_seqs[i] += comb_nts[index]
    if count_nt_combs == index_combs:
    count_nt_combs = 0
    index += 1
    if index > len(comb_nts)-1:
    index = 0

    return comb_seqs

    def number_combs(seq,stype='dna'):
    '''Counts the number of possible sequence combinations from ambiguous sequences'''

    count_combs = 1

    for nt in seq:
    if re.match(r'[M|R|W|S|Y|K]', nt, re.IGNORECASE):
    count_combs *= 2
    if re.match(r'[V|H|D|B]', nt, re.IGNORECASE):
    count_combs *= 3
    if re.match(r'[X|N]', nt, re.IGNORECASE):
    count_combs *= 4

    return count_combs

    ResponderEliminar