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)