Hola,
hoy mi colega Najla Ksouri y yo teníamos que comprobar si dos versiones del mismo genoma, una de Phytozome y otra de Ensembl Plants, eran iguales en secuencia. Para eso usamos este one-liner Perl que hace una digestión MD5 de las secuencias desde un fichero en formato FASTA, donde cada secuencia va precedida de una cabecera o header, capturada en la variable $h:
cat Prunus_persica.Prunus_persica_NCBIv2.dna.toplevel.fa | perl -MDigest::MD5=md5_hex -lne 'if(/^(>\S+)/){ $h=$1 } else { $fa{$h} .= $_ } END{ foreach $h (keys(%fa)){ print "$h\t".md5_hex(uc($fa{$h})) }}'
Éstas son las primeras líneas que imprime:
>scaffold_284 390b26eb24b1926aca0cbd2be8174139 >scaffold_112 9ed6d212a00a8018572c06d6c21fc78d >scaffold_12 4d99812ab949deb933d84a5a9ebc5787 ...
Como se puede ver, la primera columna es la primera palabra de cada cabecera, y la segunda es la digestión MD5 de esa secuencia, que sería idéntica para dos secuencias iguales.
Finalmente, veréis que llama al módulo core Digest::MD5, por tanto ya instalado, y que usa la función uc para convertir las secuencias a mayúsculas antes de digerirlas.
Hasta pronto,
Bruno
No hay comentarios:
Publicar un comentario