comprueba si dos genomas FASTA son iguales

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