mapeo de coordenadas entre ensamblajes

Hola,
una tarea con la que tropiezas a menudo cuando trabajas en genómica es la de traducir unas coordenadas de una versión del genoma a las coordenadas equivalentes de otra versión.

Ejemplo de mapeo entre dos versiones o ensamblajes del mismo genoma, tomado de http://ensemblgenomes.org/info/data/assembly_mapping.

En la figura se resume el problema. Es importante darse cuenta de que puede haber elementos en la versión n que estén ausentes en la n+1, o viceversa. Con esa salvedad, podemos traducir coordenadas en Ensembl de manera sencilla usando la interfaz REST: https://rest.ensembl.org/documentation/info/assembly_map

Por ejemplo, para convertir el fragmento 1-10 del cromosoma 1 de cebada de la versión de 2012 (ASM32608v1) a la de 2017 (IBSC_v2), podríamos hacer:

$ wget -q --header='Content-type:application/json' \

  'https://rest.ensembl.org/map/hordeum_vulgare/ASM32608v1/1:1..10/IBSC_v2?'  -O -

Obtenemos los resultados en formato JSON:

{"mappings":[{"original":{"end":10,"seq_region_name":"1","assembly":"ASM32608v1","start":1,

"strand":1,"coord_system":"chromosome"},"mapped":{"start":268488,"assembly":"IBSC_v2",

"seq_region_name":"chr1H","end":268497,"coord_system":"chromosome","strand":1}}]}

Puedes averiguar los nombres de las diferentes versiones soportadas en la herramienta Tools->AssemblyConverter (por ejemplo en humano, plantas),

hasta pronto,
Bruno

Extraer coordenadas de átomos en un fichero PDB

Los ficheros PDB (Protein Data Bank) contienen las coordenadas espaciales de los átomos de proteínas cuya estructura está resuelta por técnicas de rayos X o resonancia magnética nuclear (RMN). En una entrada de blog de Bruno podéis encontrar una muy buena explicación de estas técnicas realizada por el periódico El País.

Las estructuras de proteínas en formato PDB contienen mucha información que habitualmente no nos interesa para nuestros experimentos e interfiere con muchos programas que únicamente necesitan las coordenadas espaciales de los átomos. Además un fichero PDB normalmente contiene datos de múltiples moléculas y en diferentes posiciones en el cristal. Por ello es muy conveniente extraer de los ficheros PDB únicamente la información estructural que vamos a utilizar y descartar el resto.

En el post anterior vimos como extraer líneas de un texto o archivo. En la entrada de hoy veremos como reusar ese código para extraer los átomos deseados de un fichero PDB.

En el ejemplo se indica el fichero PDB original en la variable $pdb_file, en nuestro caso será la estructura de la famosa insulina que se puede descargar aquí (2INS.pdb). Y en el array @desired_coords le indicaremos el tipo de átomos que queremos extraer en un formato muy similar al usado por PyMOL.

   
 # File with 3D coords of insulin  
 my $pdb_file = "2INS.pdb";  
   
 # Extract coords of the alpha carbons of the chain A, and alpha and beta carbons of histidines of chain B  
 my @desired_coords = ('a/*/ca/', 'b/*/ca|cb/his');  
   
 open(PDBFILE,$pdb_file);  
 my $pdb_content = join('',<PDBFILE>);  
 close PDBFILE;  
   
 my $pdb_coords = join('',extract_pdb_coords($pdb_content,\@desired_coords));  
   
 open(PDBFILE,">$pdb_file.out");  
 print PDBFILE $pdb_coords;  
 close PDBFILE;  
   
 # Extract desired PDB coordinates from PDB entries  
 sub extract_pdb_coords {  
   
      my ($pdb_content, $types) = @_;  
   
      my $pdb_data;  
      my $patterns;  
      my @pattern_parts = ('chain','res_number','res_type','res_name');  
      foreach my $type (@{$types}) {  
           my $count = 0;  
           my %pattern;  
           while ($type =~ /([^\/]+)/g){  
                $type = $'; # Text to the right of the match  
                $pattern{$pattern_parts[$count]} = $1;  
                $count++;  
           }  
           push(@{$patterns},\%pattern);  
      }  
   
      foreach my $pattern (@{$patterns}){  
           my ($chain,$res_number,$res_type,$res_name);  
           if (!defined($pattern->{'chain'}) || !$pattern->{'chain'} || $pattern->{'chain'} eq '*'){  
                $chain = '\w{1}';  
           } else {  
                $chain = uc($pattern->{'chain'});  
           }  
           $pattern->{'res_number'} =~ s/\s+//;  
           if (!defined($pattern->{'res_number'}) || !$pattern->{'res_number'} || $pattern->{'res_number'} eq '*'){  
                $res_number = '\d+';  
           } elsif ($pattern->{'res_number'} =~ /^(\d+)$/) {  
                $res_number = $1;  
           } elsif ($pattern->{'res_number'} =~ /^(\d+)-(\d+)$/) {  
                $res_number = '('.join('|', $1 .. $2).')';  
           } elsif ($pattern->{'res_number'} =~ /\d,/) {  
                $res_number = '('.join('|', split(",",$pattern->{'res_number'})).')';  
           }  
           if (!defined($pattern->{'res_type'}) || !$pattern->{'res_type'} || $pattern->{'res_type'} eq '*'){  
                $res_type = '[\w\d]+';  
           } else {  
                $res_type = uc($pattern->{'res_type'});  
           }  
           if (!defined($pattern->{'res_name'}) || !$pattern->{'res_name'} || $pattern->{'res_name'} eq '*'){  
                $res_name = '\w{3}';  
           } else {  
                $res_name = uc($pattern->{'res_name'});  
           }  
           $pattern = '/(ATOM|HETATM)\s+\d+\s+('.$res_type.')\s+('.$res_name.')\s('.$chain.')\s+('.$res_number.')\s+.+/';  
      }  
   
      my @pdb_data_lines = extract_lines_from_text($pdb_content, $patterns);  
      if (@pdb_data_lines){  
           $pdb_data = join("\n",@pdb_data_lines);  
      }  
   
      return $pdb_data;  
   
 }  
   
 # Extract lines from text with the desired patterns  
 sub extract_lines_from_text {  
   
      my ($text, $patterns) = @_;  
   
      my @data;  
      my @lines = split("\n",$text);  
   
      foreach my $line (@lines){  
           foreach my $pattern (@{$patterns}){  
                if ($pattern =~ /^\/(.+)\/$/){  
                     if ($line =~ /$1/){  
                          push(@data,$line);  
                          last;  
                     }  
                } else {  
                     if ($line =~ /\Q$pattern\E/){  
                          push(@data,$line);  
                          last;  
                     }  
                }  
           }  
      }  
   
      return @data;  
   
 }