28 de agosto de 2017

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

1 comentario:

  1. La expresión regular se puede quitar, ya que se trata de localizar solo un carácter en una determinada posición.

    Se puede reducir a simplemente /^>/, pero hay otra opción ligeramente más rápida: if (">" eq substr $_, 0, 1)

    ResponderEliminar