Mostrando entradas con la etiqueta bioperl. Mostrar todas las entradas
Mostrando entradas con la etiqueta bioperl. Mostrar todas las entradas

2 de junio de 2022

Traducir codones a aminoácidos en el terminal

Hola,

hoy me he vuelto a encontrar con una tarea frecuente cuando trabajas con secuencias genómicas, la de traducir codones a aminoácidos. Para unas pocas secuencias se pueden usar herramientas Web como https://web.expasy.org/translate , pero para hacerlo de manera programática y a medida me pareció más conveniente hacerlo con el clásico módulo BioPerl, que tiene documentación específica para esta tarea aquí, incluyendo cómo cambiar de tabla de uso de codones o de fase. En concreto usaremos un one-liner, un pequeño comando en Perl, similar a otros que podeís ver aquí.

 



##1) Instalación, muestro 3 opciones


# 1.1) ubuntu/Debian
sudo apt-get install -y bioperl

# 1.2) bioconda
conda install perl-bioperl

# 1.3) cpanm
cpan -i App:cpanminus cpanm --force Bio::Perl


# 2) Traducciones de un fichero FASTA (codons.fna)

# 2.1 Si bioperl está en ubicación central en @INC
perl -MBio::Seq -lne 'if(/^>/){ print } else { print Bio::Seq->new(-seq => $_)->translate()->seq() }' codons.fna > peptides.faa 


# 2.2) Si bioperl está instalada en otro sitio
perl -I $biopath -MBio::Seq -lne 'if(/^>/){ print } else { print Bio::Seq->new(-seq => $_)->translate()->seq() }' codons.fna > peptides.faa

# 2.3) Para comprobar codones de stop internos
grep -B 1 -P "\*[A-Z]" peptides.faa


Hasta pronto,

Bruno

4 de julio de 2017

rooting and laddering Newick trees

Buenas,
esta entrada es para compartir un script Perl para enraizar árboles, muchos árboles, de manera automática, y ordenar los nodos por distancia. Pongo el título en inglés para los buscadores.
Rubén Sancho y yo estamos probando el software GRAMPA, que requiere una colección de árboles de genes precalculados en formato Newick, pero además los quiere enraizados y ordenados de manera ascendente. Como son más de mil árboles no lo queríamos hacer uno a uno con FigTree, por poner un ejemplo; queríamos automatizarlo y para ello aprovechamos los módulos Bio::TreeIO y Bio::Phylo. El primero es de Bioperl y ya lo tenía instalado, el segundo lo instalé con: $ sudo cpanm Bio::Phylo

El árbol de ejemplo, contenido en el fichero unrooted.ph, es:

(1_Sbic:0.047,2_Osat:0.068,(((((((((((3_B422:0.007,(((4_Barb:0.000,((5_Bret:0.003,(6_BsyE:0.000,7_BsyE:0.000):0.000):0.000,8_BsyC:0.000):0.000):0.000,9_Bpho:0.002):0.000,(10_BsyC:0.000,11_BsyC:0.000):0.001):0.000):0.008,((12_B422:0.005,13_Bpho:0.005):0.002,14_Barb:0.000):0.005):0.000,((((15_Bdis:0.000,16_Bdis:0.000):0.000,17_Bmex:0.033):0.013,18_Bhyb:0.001):0.014,(19_Bboi:0.021,((20_Bmex:0.017,21_Bsta:0.034):0.012,22_Bsta:0.000):0.003):0.020):0.012):0.002,23_Barb:0.000):0.000,24_Brup:0.012):0.000,25_BsyG:0.000):0.000,(26_Bpin:0.000,27_BsyG:0.000):0.000):0.011,28_Bsta:0.000):0.053,29_BsyG:0.000):0.036,(30_Barb:0.005,(31_Bpho:0.024,(32_BsyC:0.000,33_BsyE:0.000):0.027):0.029):0.000):0.005,34_Bboi:0.030):0.035);

Invocando el script lo enraizamos con el taxón elegido, el outgroup '_Sbic':

$ perl reroot_tree.pl unrooted.ph > rooted.ph

Obtenemos el fichero rooted.ph, que podemos visualizar con FigTree:

Árbol enraizado y con orden creciente de nodos.
Éste es el código de reroot_tree.pl :

 #!/usr/bin/perl -w  
   
 # Re-roots an input Newick tree with a user-defined outgroup and   
 # prints the resulting tree in ascending or descending node order.  
 # Based on https://github.com/phac-nml/snvphyl-tools/blob/master/rearrange_snv_matrix.pl  
   
 use strict;  
 use Bio::TreeIO;  
 use Bio::Phylo::IO;  
 use Bio::Phylo::Forest::Tree;  
    
 my $OUTGROUPSTRING = '_Sbic'; # change as needed  
 my $NODEORDER = 1; # 1:increasing, 0:decreasing  
 my ($outfound,$outnode,$sorted_newick) = (0);  
   
 die "# usage: $0 <tree.newick>\n" if(!$ARGV[0]);  
   
 # read input tree  
 my $input = new Bio::TreeIO(-file=>$ARGV[0],-format=>'newick');  
 my $intree = $input->next_tree();  
   
 # find outgroup taxon  
 for my $node ( $intree->get_nodes() ) {   
  if(defined($node->id()) && $node->id() =~ m/$OUTGROUPSTRING/) {  
   $outnode = $node;   
   $outfound = 1;  
   last;  
  }  
 }  
 if($outfound == 0) {  
  die "# cannot find outgroup $OUTGROUPSTRING in input tree $ARGV[0]\n";  
 }  
   
 # root in outgroup and sort in increasing order  
 $intree->reroot_at_midpoint($outnode);  
   
 # sort nodes in defined order and print  
 my $unsorted_tree = Bio::Phylo::IO->parse(  
  '-string' => $intree->as_text('newick'),  
  '-format' => 'newick'  
 )->first();  
   
 $unsorted_tree->ladderize($NODEORDER);  
   
 $sorted_newick = $unsorted_tree->to_newick();  
 $sorted_newick =~ s/'//g;  
   
 print $sorted_newick;  


Un saludo,
Bruno

3 de noviembre de 2015

Anotando datos del NCBI de forma automática con un script de Perl

Esta semana me di cuenta que un antiguo script de Perl que leía datos bibliográficos de Pubmed automáticamente había dejado de funcionar.

La explicación es que el NCBI ha dejado de dar soporte a su servicio web con el protocolo SOAP desde el 1 de julio de 2015. No soy un experto informático, pero hace unos años ese protocolo estaba de moda y ahora parece que ya no (imagino que por razones de seguridad).

He aquí el código OBSOLETO que da el error:

 use SOAP::Lite;  
   
 my $pubmed_id = '24234003';  
   
 my $WSDL = 'http://www.ncbi.nlm.nih.gov/soap/v2.0/efetch_pubmed.wsdl';  
 my $soap = SOAP::Lite->service($WSDL)  
                ->on_fault( sub { my($self, $res) = @_;  
                     die  
                     "faultcode:", ref $res ? $res->faultcode : $self->transport->status, "\n" ,  
                     "faultstring:", ref $res ? $res->faultstring : $self->transport->status, "\n";  
                });  
 my $response = $soap->run_eFetch(     SOAP::Data->name(id => $pubmed_id),  
                          SOAP::Data->name(db => "pubmed")  
                     );  
   
 my %results;  
   
 if (ref($response) eq "HASH") {  
      $results{'title'} = $response->{'PubmedArticle'}{'MedlineCitation'}{'Article'}{'ArticleTitle'};  
      $results{'journal'} = $response->{'PubmedArticle'}{'MedlineCitation'}{'Article'}{'Journal'}{'Title'};  
      $results{'pubmed'} = $response->{'PubmedArticle'}{'MedlineCitation'}{'PMID'};  
      $results{'volume'} = $response->{'PubmedArticle'}{'MedlineCitation'}{'Article'}{'Journal'}{'JournalIssue'}{'Volume'};  
      $results{'issue'} = $response->{'PubmedArticle'}{'MedlineCitation'}{'Article'}{'Journal'}{'JournalIssue'}{'Issue'};  
      $results{'year'} = $response->{'PubmedArticle'}{'MedlineCitation'}{'Article'}{'Journal'}{'JournalIssue'}{'PubDate'}{'Year'};  
      $results{'pages'} = $response->{'PubmedArticle'}{'MedlineCitation'}{'Article'}{'Pagination'}{'MedlinePgn'};  
      $results{'url'} = 'http://www.ncbi.nlm.nih.gov/pubmed/'.$results{'pubmed'};  
      if (ref($response->{'PubmedArticle'}{'MedlineCitation'}{'Article'}{'AuthorList'}{'Author'}) eq "ARRAY") {  
           my @authors;  
           foreach my $author (@{$response->{'PubmedArticle'}{'MedlineCitation'}{'Article'}{'AuthorList'}{'Author'}}){  
                my $initials = join('.',split('',$author->{'Initials'}));  
                push(@authors, $author->{'LastName'}." ".$initials);  
           }  
           $results{'authors'} = join(", ",@authors);  
      } elsif (ref($response->{'PubmedArticle'}{'MedlineCitation'}{'Article'}{'AuthorList'}{'Author'}) eq "HASH") {  
           my $author = $response->{'PubmedArticle'}{'MedlineCitation'}{'Article'}{'AuthorList'}{'Author'};  
           my $initials = join('.',split('',$author->{'Initials'}));  
           $results{'authors'} = $author->{'LastName'}." ".$initials;  
      }  
      foreach my $key (keys %{$results}){  
           if (!defined($results{$key})) {  
                $results{$key} = '';  
           }  
      }  
 }  
   
 print %results;  
   


Ahora las consultas automáticas al NCBI se han de realizar insertando nuestra consulta en una URL e indicando el formato deseado de respuesta. Por ejemplo, la consulta equivalente al código SOAP anterior se puede realizar con la siguiente URL:
 O si preferimos el resultado en XML:
De esta forma, el código anterior queda reescrito de la siguiente forma:

 use LWP::Simple;  
   
 my $pubmed_id = '24234003';  
   
 my $eutils_url = 'http://eutils.ncbi.nlm.nih.gov/entrez/eutils/';  
 my $eutil = 'efetch';  
 my $db = 'pubmed';  
 my $type = 'docsum';  
 my $mode = 'text';  
   
 my $url = sprintf('%s%s.fcgi?db=%s&id=%s&rettype=%s&retmode=%s', $eutils_url, $eutil, $db, $pubmed_id, $type, $mode);  
   
 my $data = get($url);  
   
 my %results;  
   
 if (defined($data)){  
      $data =~ s/\012\015?|\015\012?/ /g;  
      if ($data =~ /\d+\:\s*(.+?)\.(.+?)\.(.+?)\.(.+?)\s.+?\;(\d+)\((\d+)\)\:([\d\-]+)./){  
           @{$results{'authors'}} = split(/,\s?/,$1);  
           $results{'title'} = $2;  
           $results{'journal'} = $3;  
           $results{'year'} = $4;  
           $results{'volume'} = $5;  
           $results{'issue'} = $6;  
           $results{'pages'} = $7;  
           $results{'pubmed'} = $pubmed_id;  
           $results{'url'} = 'http://www.ncbi.nlm.nih.gov/pubmed/'.$pubmed_id;  
      } elsif ($data =~ /\d+\:\s*(.+?)\.(.+?)\.(.+?)\.(.+?)\s/){  
           @{$results{'authors'}} = split(/,\s?/,$1);  
           $results{'title'} = $2;  
           $results{'journal'} = $3;  
           $results{'year'} = $4;  
           $results{'pubmed'} = $pubmed_id;  
           $results{'url'} = 'http://www.ncbi.nlm.nih.gov/pubmed/'.$pubmed_id;  
      }  
 }  
   
 while (my ($key, $value) = each %results) {  
      if (ref($value) ne ARRAY){  
           printf("%s: %s\n", $key, $value);  
      } else {  
           printf("%s: %s\n", $key, join(', ',@$value));  
      }  
 }  
   

Para otro tipo de consultas automáticas (secuencias, artículos científicos...) recomiendo leer la documentación oficial del NCBI, donde se explica la utilidad de las diferentes Entrez Programming Utilities (E-utilities): ESearch, EFetch, ESummary, EPost y ELink.

Seguro que alguno de vosotros está pensando porqué no usar BioPerl o BioPython para estas consultas. La respuesta es porque quería explicar el fin de la era SOAP y la estructura de las nuevas consultas vía URL.

Aquí el código en BioPerl:

 use Bio::DB::EUtilities;  
   
 my $pubmed_id = '24234003';  
   
 my $eutil = 'efetch';  
 my $db = 'pubmed';  
 my $type = 'docsum';  
 my $mode = 'text';  
 my $email = 'mymail@foo.bar';  
    
 my $factory = Bio::DB::EUtilities->new(     -eutil  => $eutil,  
                          -db   => $db,  
                          -id   => $pubmed_id,  
                          -email  => $email,  
                          -rettype => $type,  
                          -retmode => $mode  
                     );  
    
 my $data = $factory->get_Response->content;  
   
 my %results;  
   
 if (defined($data)){  
      $data =~ s/\012\015?|\015\012?/ /g;  
      if ($data =~ /\d+\:\s*(.+?)\.(.+?)\.(.+?)\.(.+?)\s.+?\;(\d+)\((\d+)\)\:([\d\-]+)./){  
           @{$results{'authors'}} = split(/,\s?/,$1);  
           $results{'title'} = $2;  
           $results{'journal'} = $3;  
           $results{'year'} = $4;  
           $results{'volume'} = $5;  
           $results{'issue'} = $6;  
           $results{'pages'} = $7;  
           $results{'pubmed'} = $pubmed_id;  
           $results{'url'} = 'http://www.ncbi.nlm.nih.gov/pubmed/'.$pubmed_id;  
      } elsif ($data =~ /\d+\:\s*(.+?)\.(.+?)\.(.+?)\.(.+?)\s/){  
           @{$results{'authors'}} = split(/,\s?/,$1);  
           $results{'title'} = $2;  
           $results{'journal'} = $3;  
           $results{'year'} = $4;  
           $results{'pubmed'} = $pubmed_id;  
           $results{'url'} = 'http://www.ncbi.nlm.nih.gov/pubmed/'.$pubmed_id;  
      }  
 }  
   
 while (my ($key, $value) = each %results) {  
      if (ref($value) ne ARRAY){  
           printf("%s: %s\n", $key, $value);  
      } else {  
           printf("%s: %s\n", $key, join(', ',@$value));  
      }  
 }  
   

Y el código en BioPython:

 from Bio import Entrez  
   
 pubmed_id = '24234003';  
   
 handle = Entrez.efetch(db='pubmed', id=pubmed_id, rettype='docsum', retmode="xml", email='mymail@foo.bar')  
   
 results = Entrez.read(handle)  
   
 for key, value in results[0].iteritems():  
      print "%s: %s" % (key,value)  
   



11 de enero de 2013

Indexando archivos FASTA

Hola y feliz año nuevo con retraso!
Recientemente me he visto en la necesidad de acceder de maneara aleatoria a secuencias de un archivo FASTA de gran tamaño, en mi caso cercano a los 2 Gb,
con las secuencias formateadas en una sola línea.

Primero probé con las herramientas de BLAST+, en concreto:
$ makeblastdb -in archivo.fasta -parse_seqid

Con la idea de posteriormente consultar las secuencias con algo como (más ejemplos aquí):

$ blastdbcmd -query -db archivo.fasta

Sin embargo, la generación del índice se hizo eterna y de hecho la aborté, por problemas que desconozco y no he seguido mirando.

Posteriormente, tras buscar en Google encontré dos caminos en Perl:

1)  con ayuda de scripts de Bioperl: bp_index.pl y bp_fetch.pl

2) con ayuda del módulo core Tie::File, como muestro en el ejemplo a continuación, mi solución preferida, un saludo:

 #!/usr/bin/perl -w  
 use strict;  
 use Tie::File;  
   
 my $fasta_file = '/path/to/file.fasta';  
   
 my (%index_fasta,@fasta_array);  
 open(FASTA,$fasta_file) ||   
    die "# cannot open $fasta_file\n";  
 while(<FASTA>)  
 {  
    if(/^> any typical marker in header, such as gi (\S+)/)  
    {  
       $index_fasta{$1} = $.;   
    }  
 }  
 close(FASTA);   
 printf("# indexed %d sequences in file %s\n\n",  
    scalar(keys(%index_fasta)),$fasta_file);  
   
 tie(@fasta_array,'Tie::File',$fasta_file) ||   
    die "# cannot tie file $fasta_file\n";  
   
 print "ejemplo: secuencia con identificador 'id12':\n";  
 print "$fasta_array[$index_fasta{'id12'}-1]\n".  
    "$fasta_array[$index_fasta{'id12'}]\n";  

1 de junio de 2010

Extraer secuencias intergénicas de archivos GenBank

La tarea de hoy consiste en leer un archivo GenBank, por ejemplo de un genoma circular bacteriano, con el fin de extraer las secuencias intergénicas. Para esta tarea vamos a utilizar código de BioPerl y antes necesitaremos descargar al menos un fichero GenBank, como se explica por ejemplo en el taller de (bio)perl. Una vez tengamos los archivos de entrada necesarios, podemos resolver el problema de las secuencias intergénicas con una subrutina como ésta:

1:  use Bio::SeqIO;
2:
3: # prog1.pl
4: # recibe hasta 4 argumentos:
5: # 1) archivo entrada en formato GenBank
6: # 2) nombre del archivo de salida en formato FASTA
7: # 3) longitud mínima de las regiones intergénicas (opcional, por defecto toma todas)
8: # 4) longitud máxima de las regiones intergénicas (opcional)
9:
10: sub extract_intergenic_from_genbank
11: {
12: my ($infile,$out_intergenic_file,$min_intergenic_size,$max_intergenic_size) = @_;
13:
14: my ($n_of_intergenic,$gi,$start,$end,$length,$strand,$taxon) = (0);
15: my $in = new Bio::SeqIO(-file => $infile, -format => 'genbank' );
16:
17: open(FNA,">$out_intergenic_file") ||
18: die "# extract_intergenic_from_genbank : cannot create $out_intergenic_file\n";
19:
20: while( my $seq = $in->next_seq) # scan all sequences found in $input
21: {
22: my ($gbaccession,$sequence,$gen,@genes) = ( $seq->accession() );
23: $sequence = $seq->primary_seq()->seq() ||
24: 'no encuentro la secuencia primaria, necesito un fichero GenBank completo!';
25: $taxon = '';
26:
27: for my $f ($seq->get_SeqFeatures)
28: {
29: if($f->primary_tag =~ /CDS|rRNA|tRNA/) # campos de 'genes'
30: {
31: $gi = '';
32: if($f->has_tag('db_xref'))
33: {
34: my $crossrefs = join(',',sort $f->each_tag_value('db_xref'));
35: if($crossrefs =~ /(GI\:\d+)/){ $gi = $1 }
36: }
37: elsif($f->has_tag('locus_tag'))
38: {
39: if($gi eq '' && $f->has_tag('locus_tag'))
40: {
41: $gi = "ID:".join(',',sort $f->each_tag_value('locus_tag'));
42: }
43: }
44:
45: next if($gi eq '');
46:
47: $start = $f->start();
48: $end = $f->end();
49: $strand = $f->location()->strand();
50: push(@genes,[$start,$end,$gi,$strand]);
51: }
52: elsif($f->primary_tag() =~ /source/)
53: {
54: if($f->has_tag('organism'))
55: {
56: foreach my $element ($f->each_tag_value('organism'))
57: {
58: $taxon .= "[$element],";
59: }
60: chop $taxon;
61: }
62: }
63: }
64:
65: # calcular ORFs vecinos y corta las secuencias intergenicas
66: for($gen=1;$gen<scalar(@genes);$gen++)
67: {
68: next if($genes[$gen-1][1] >= $genes[$gen][0]); #no solapa asume genes en orden
69: next if( defined($min_intergenic_size) && $min_intergenic_size > 0 &&
70: $genes[$gen][0]-$genes[$gen-1][1]+1 < $min_intergenic_size); # evita cortas
71: next if( defined($max_intergenic_size) && $max_intergenic_size > 0 &&
72: $genes[$gen][0]-$genes[$gen-1][1]+1 > $max_intergenic_size); # evita largas
73:
74: $n_of_intergenic++;
75:
76: # calcula bordes de intergenes
77: $start = $genes[$gen-1][1] + 1;
78: $end = $genes[$gen][0] - 1;
79: $length = $end - $start + 1;
80:
81: print FNA ">intergenic$n_of_intergenic|$gbaccession|coords:$start..$end|".
82: "length:$length|neighbours:$genes[$gen-1][2]($genes[$gen-1][3]),".
83: "$genes[$gen][2]($genes[$gen][3])|$taxon\n";
84: print FNA substr($sequence,$start-1,$length)."\n";
85: }
86: }
87:
88: close(FNA);
89:
90: return $n_of_intergenic;
91: }