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
La expresión regular se puede quitar, ya que se trata de localizar solo un carácter en una determinada posición.
ResponderEliminarSe puede reducir a simplemente /^>/, pero hay otra opción ligeramente más rápida: if (">" eq substr $_, 0, 1)