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

contenido GC de un fichero FASTA

Hola,

una pregunta habitual cuando analizas un fichero de nucleótidos, por ejemplo un ensamblaje de un genoma, es qué porcentaje GC tiene.  Asumiendo que el fichero está en formato FASTA, podemos obtener fácilmente ese valor con un mini-programa (one-liner) escrito en lenguaje perl. Por ejemplo, para el genoma comprimido de Brachypodium distachyon obtenido de Ensembl Plants, podríamos obtenerlo así:

zcat Brachypodium_distachyon.Brachypodium_distachyon_v3.0.dna.toplevel.fa.gz | \
   perl -lne 'if(!/^>/){ $SQ=uc($_); while($SQ =~ /([ACTG])/g){ $stat{$1}++; $tot++ } } 
   END{ printf("%%GC=%1.1f\n",100*($stat{"G"}+$stat{"C"})/$tot);  
      foreach $nt (keys(%stat)){ print "$nt\t$stat{$nt}" } }'

%GC=46.4
A	72549289
T	72561114
C	62839311
G	62789747

Si quieres calcular el %GC solamente para ciertas regiones del genoma entonces puedes codificarlas en un fichero BED y usar bedtools nuc, como se explica en https://www.biostars.org/p/47047

Hasta pronto,

Bruno

más one-liners Perl

Hola,
antes de que se acabe el año aprovecho para compartir con vosotros un excelente tutorial de one-liners de Perl, esos comandos que en una línea permiten ejecutar complejas operaciones en el terminal de Linux, el símbolo del sistema de Windows, o, mejor aún, desde dentro de una ventana de MobaXterm.
El tutorial se aloja e:

https://github.com/learnbyexample/Command-line-text-processing/blob/master/perl_the_swiss_knife.md

y tiene ejemplos tan útiles como:

# 1) calcula máximo de una lista de números separados por comas

$ echo '34,17,6' | perl -MList::Util=max -F, -lane 'print max @F'
34

# 2) valida y expande un one-liner a un programa completo más comprensible

$perl -MO=Deparse -ne 'if(!$#ARGV){$h{$_}=1; next} print if $h{$_}'
LINE: while (defined($_ = )) {
    unless ($#ARGV) {
        $h{$_} = 1;
        next;
    }
    print $_ if $h{$_};
}
-e syntax OK

El tutorial tiene también recetas para usar el resto de herramientas del terminal Linux, como grep, sed y muchas otras, en

https://github.com/learnbyexample/Command-line-text-processing

Feliz año!
Bruno

Secuencia de referencia para experimento TagSeq

Hola,
cada vez se van publicando más trabajos donde se emplea TagSeq, una versión low cost de RNAseq que se especializa en secuenciar el máximo número de transcritos posibles, pero sólo unos cuantos cientos de bases de su extremo 3', contando desde la cola poliA. Un tamaño típico de librería TagSeq es 500b.

Protocolo TagSeq, tomado de https://tinyurl.com/y9yc4u5a.

Cuando obtenemos lecturas o reads de este tipo y las queremos alinear contra los transcritos anotados del genoma de referencia puede ser útil, con vistas a posibles normalizaciones posteriores que consideren la longitud original del gen, recortar las secuencias de referencia. Os pongo un ejemplo en Perl:

zcat primaryTranscriptOnly.fa.gz | \
     perl -lne 'if(/^(>.*)/){$h=$1}else{$fa{$h} .= $_} END{ foreach $s (sort keys(%fa)){ print "$s\n".substr($fa{$s},-500)."\n" }}' > \     
     primaryTranscriptOnly.TagSeq500b.fa

Hasta luego,
Bruno

mean sequence length in FASTA

Hola,
para calcular la longitud promedio de las secuencias de un fichero FASTA estábamos usando el siguiente comando de Perl en el terminal, modificado de otro en https://eead-csic-compbio.github.io/perl_bioinformatica/node90.html:

$ time perl -lne 'if(/^(>.*)/){$h=$1}else{$fa{$h}.=$_} END{ foreach $h (keys(%fa)){$m+=length($fa{$h})}; printf("%1.0f\n",$m/scalar(keys(%fa))) }' ~/db/swissprot
376

real 0m4.943s
user 0m4.488s
sys 0m0.452s

Sin embargo, para ficheros de metagenomas era muy lento, y lo hemos reemplazado por este otro, que no guarda las secuencias en memoria (en un hash llamado %fa):

$ time perl -lne 'if(/^(>.*)/){$nseq++;$m+=$l;$l=0}else{$l+=length($_)} END{ $m+=$l; printf("%1.0f\n",$m/$nseq) }' ~/db/swissprot

376

real 0m2.013s
user 0m1.880s
sys 0m0.132s

Hasta luego,
Bruno

Actualización de Perl en bioinformática

Buenas,
solamente quería comentar que el ya vetusto curso de Perl en bioinformático ha cambiado un poco de cara y ha sido revisado para traerlo al 2017, con algunas referencias a Perl6 y nuevos one-liners, que es lo que usaremos de este material en el máster de biotecnología cuantitativa de la Universidad de Zaragoza a partir del curso que viene.

Lo podéis encontrar en: https://eead-csic-compbio.github.io/perl_bioinformatica

Portada de un libro de one-liners, tomada de http://www.catonmat.net/blog/perl-one-liners-explained-part-one

Buen finde,
Bruno

One-liner para calcular el N50 de un ensamblaje

Hola,
el estadístico N50 se usa a menudo para resumir a grosso modo un ensamblaje de secuencias, que generalmente es un fichero FASTA con una serie de contigs.
Modificando la definición de la wikipedia, N50 es la longitud de contigs tal que usando contigs de igual o mayor tamaño recuperamos la mitad de las bases de ese ensamblaje.

Figura tomada de http://bmpvieira.github.io/assembly14.

Hoy solamente quería compartir un comando one-liner de Perl que nos permite calcularlo:

$ perl -lne 'if(/^>(\S+)/){$h=$1} else {$l=length($_);$TL+=$l;$L{$h}+=$l} \

END{ foreach $s (sort {$L{$b}<=>$L{$a}} keys(%L)){ $t+=$L{$s}; \

if($t>$TL/2){ print "N50=$L{$s}";exit }}}' assembly.fasta

Cuidado al copiarlo al terminal desde el blog, mejor hacerlo línea a línea,
hasta pronto, Bruno

clustal one-liner con parallel

Hola,
hoy comparto un comando que a veces utilizo cuando necesito calcular muchos alineamientos múltiples a partir de una colección de ficheros de secuencias en formato FASTA. Como mi máquina, igual que la de casi todos, tiene amplia RAM y muchos cores, es un trabajo ideal para parallel. Supongamos que los archivos de salida están en la carpeta 'entrada' y queremos guardar los ficheros de salida en la carpeta 'path/to/salida', y que tenemos 20 cores disponibles:

$ mkdir /path/to/salida/

$ cd entrada

$ ls -1 *fasta | parallel --gnu -j 20 ~/soft/clustal-omega-1.2.1/src/clustalo \

--threads=1 -i {} -o /path/to/salida/{} :::

Este comando pondrá a trabajar 20 cores del sistema hasta que todos los archivos FASTA de la carpeta entrada estén alineados, con ganancias de tiempo de ejecución importantes en un experimento con 100 ficheros:

| cores (-j) | time(real) | time(user) | time(sys) |
|   1        | 4m34.440s  | 4m5.180s   | 0m2.168s  |
|  10        | 0m29.358s  | 3m57.768s  | 0m2.400s  |
|  20        | 0m23.248s  | 5m6.204s   | 0m3.364s  |

Un saludo,
Bruno

R one-liners

Hola,
en esta entrada quería compartir unos ejemplos para aprovechar las capacidades del lenguaje R para analizar de manera rápida datos desde el terminal, espoleado por el primero de ellos, que me pasó Carlos Cantalapiedra la semana pasada. Pero vayamos al grano usando como ejemplo un fichero de salida tabular de BLAST, al que llamarmos 'all.blast'.

Receta 1. Tienes un fichero con columnas de números y quieres saber la media de una de ellas, por ejemplo la tercera. Para la la desviación estándar cambia 'mean' por 'sd':

$ cut -f 3 all.blast | Rscript -e 'data=scan(file="stdin"); mean(data)'

Receta 2. Tienes un fichero con columnas de números y quieres calcular un histograma de una de ellas, por ejemplo la tercera:

$ cut -f 3 all.blast | \

Rscript -e 'data=scan(file="stdin"); pdf("hist.pdf");hist(data)' 

Se puede modificar mostrando el PDF creado inmediatamente:

$ cut -f 3 all.blast | \

Rscript -e 'data=scan(file="stdin"); pdf("hist.pdf");hist(data)'; evince hist.pdf 

Receta 3. Tienes un fichero con columnas de números y quieres calcular boxplots con dos de ellas, por ejemplo la primera y la segunda:

cut -f 1,2 all.blast | \

Rscript -e 'd=matrix(scan(file="stdin"),ncol=2); pdf("box.pdf");boxplot(d[,1],d[,2])'

Espero que os sirvan de inspiración,
un saludo,
Bruno