#!/usr/bin/perl
#
# Ortólogos
#
# Joaquín Ferrero, 20140406
#
# Dada una serie de archivos que contienen una relación de genes de organismos,
# y otra serie de archivos que indican los ortólogos entre los genes de distintos
# organismos, obtener una tabla con el resumen de ortólogos entre todos los organismos.
#
use v5.14;
use autodie;
use File::Basename;
use File::Slurp;
### Constantes
my $DIR = './ortologos';
### Leer archivos con datos ####################################################
my @nombres_archivos = glob("$DIR/*.id");
my @organismos; # nombres de los organismos
my %genes_de_organismo; # genes por cada organismo
for my $archivo (@nombres_archivos) {
my $organismo = basename($archivo, '.id');
push @organismos, $organismo; # el nombre es a partir del archivo
my @genes = read_file($archivo, chomp => 1); # leemos sus genes
$genes_de_organismo{$organismo} = [ @genes ]; # y los guardamos
}
@organismos = sort @organismos; # ordenación alfabética
# traductor, de nombre de organismo a número de columna dentro de la tabla
my %organismo_a_columna = map { $organismos[$_] => $_ } 0 .. $#organismos;
### Leer archivos de relaciones ################################################
my %ortologos;
for my $archivo (glob("$DIR/*.comp")) { # para todos los archivos
my($org1, $org2) = $archivo =~ /(\w+)_vs_(\w+)/; # organismos que se relacionan
my @relaciones = read_file($archivo, chomp => 1); # leemos relaciones
for (@relaciones) {
my($gen1, $gen2) = split; # los dos genes que intervienen
## La estructura a crear es la de un hash de arrays
##
## %ortologos{ ortologo } => [ ortologo [, ortologo, ...] ]
##
## Por cada ortólogo, guardamos un array anónimo con el resto de ortólogos
## con los que se relaciona.
##
## Los ortólogos tienen la forma: "$organismo;$gen"
##
push @{$ortologos{"$org1;$gen1"}}, "$org2;$gen2"; # creamos dos relaciones:
push @{$ortologos{"$org2;$gen2"}}, "$org1;$gen1"; # (relación conmutativa)
}
}
### Análisis ###################################################################
my %visto; # bandera para evitar repeticiones
for my $organismo (@organismos) { # para todos los organismos
for my $gen (@{$genes_de_organismo{$organismo}}) { # para todos sus genes
next if $visto{"$organismo;$gen"}++; # excepto si ya lo hemos visto
# (caso de un gen que ya participaba en un ortólogo)
## Iniciamos un bucle que pintará todas las relaciones encontradas
## para ese organismo y gen. Debemos repetirlo mientras existan ortólogos posibles para ese gen.
## Es importante que sea do(), y no while(): es para los casos de genes sueltos (sin relación con los demás,
## el bucle debe ejecutarse al menos una vez, para que salga en salida).
do {
## Buscamos la mayor relación de ortólogos, a partir de ese gen, de ese organismo
my @ortologos_encontrados = busca_maximo_ortologos("$organismo;$gen");
## Pintar los ortólogos encontrados
my @fila = ('') x @organismos; # inicializamos a columnas vacías
for my $ort (@ortologos_encontrados) { # para todos los encontrados
my($org,$gen) = split /[;]/, $ort;
$fila[$organismo_a_columna{$org}] = $gen; # ponemos el gen en la columna correspondiente
}
say join '', map { sprintf "%-11s", $_ } @fila; # y pintamos
## Borrar los ortólogos encontrados.
## Debemos destruir los ortólogos que hemos encontrado y pintado,
## para que no influyen en el resto de búsquedas (proceso destructivo)
## El proceso consiste en eliminar los ortólogos que hemos encontrado
## de los array de ortólogos relacionados con cada ortólogo encontrado
## (esto no queda muy claro, ¿verdad?)
my $filtro = join '|', map{ quotemeta } @ortologos_encontrados; # patrón de exp. reg.
for my $ortologo (@ortologos_encontrados) { # para todos los ortólogos encontrados
my($org,$gen) = split /[;]/, $ortologo;
next if not exists $ortologos{$ortologo}; # caso especial:
# ese ortólogo ya fue borrado antes
# le quitamos a ese ortólogo, todos los ortólogos con los que se relaciona,
# que intervienen en el conjunto de ortólogos encontrados
$ortologos{$ortologo} = [ grep {! /$filtro/} @{$ortologos{$ortologo}} ];
# Si ese ortólogo ya no se relaciona con más ortólogos
if (@{$ortologos{$ortologo}} == 0) {
delete $ortologos{$ortologo}; # lo borramos
$visto{$ortologo}++; # lo marcamos para no procesarlo más
}
}
} while ($ortologos{"$organismo;$gen"}); # y repetimos hasta que no haya más relaciones
# de ese ortólogo
}
}
### Proceso recursivo de búsqueda de la máxima cantidad de ortólogos relacionados
## Vamos llamando a la propia función con un ortólogo más cada vez, hasta que no puede más.
## Como argumento recibe un array con los ortólogos relacionados encontrados hasta ahora.
## Devuelve la mayor relación encontrada, o la misma que recibió como argumento, si no puede mejorarla.
sub busca_maximo_ortologos {
my @ortologos_temporal = @_;
## Paso 1. Ver si el último ortólogo temporal tiene relación con todos los ortólogos anteriores
my $ultimo_ortologo = $ortologos_temporal[-1];
if (not exists $ortologos{$ultimo_ortologo}) { # caso de no existir ninguna relación
return @ortologos_temporal; # regresamos con lo que tenemos
}
my $nrelaciones = 0; # contaremos el número de relaciones entre ortólogos
my %organismos_vistos; # cada organismo contará una sola vez
# sacamos el listado de ortólogos con los que se relaciona el último ortólogo
my %ortologos_ultimo_ortologo = map { $_ => 1 } @{$ortologos{$ultimo_ortologo}};
for my $i (0 .. $#ortologos_temporal-1) { # para todos los ortólogos anteriores al último
my $ortologo_anterior = $ortologos_temporal[$i];
my($org,$gen) = split /[;]/, $ortologo_anterior;
# si es de un organismo que ya cuenta con un ortólogo entre los ortólogos que estamos analizando
if ($organismos_vistos{$org}++) {
delete $ortologos_ultimo_ortologo{$ortologo_anterior}; # no lo vamos a visitar
}
# si existe una relación entre el último ortólogo, y el anterior
if ($ortologos_ultimo_ortologo{$ortologo_anterior}) {
$nrelaciones++; # contamos una relación mas
delete $ortologos_ultimo_ortologo{$ortologo_anterior}; # y no lo visitaremos (ya lo fue)
}
}
# para que el último ortólogo sea legal, debe tener tantas relaciones como ortólogos anteriores a él
if ($nrelaciones != @ortologos_temporal-1) { # no hay relaciones completas con los anteriores
return @ortologos_temporal[0 .. $#ortologos_temporal-1]; # devolvemos todos menos el último, por ser feo
}
## Paso 2. Probar con los ortólogos relacionados con el último ortólogo
# Hacemos un bucle por todos los ortólogos con los que se relaciona el último ortólogo,
# y que todavía no han sido visitados
my @maximo_ortologos; # Variables para recordar el máximo encontrado
for my $prueba_ortologo (keys %ortologos_ultimo_ortologo) { # para todos los que quedan por visitar
# vemos si ese ortólogo de prueba marca un nuevo record
my @ortologos_encontrados = busca_maximo_ortologos(@ortologos_temporal, $prueba_ortologo);
if (@maximo_ortologos < @ortologos_encontrados) { # sí
@maximo_ortologos = @ortologos_encontrados; # lo recordaremos
last if @maximo_ortologos == @organismos; # caso de encontrar un ortólogo máximo (extremo)
# no hace falta seguir buscando
}
}
# Paso 3. Elegir la combinación mayor de lo encontrado hasta ahora
if (@maximo_ortologos) { # si hemos superado lo encontrado hasta ahora
return @maximo_ortologos; # devolvemos el récord de ortólogos encontrados
}
else { # si no
return @ortologos_temporal; # devolvemos lo encontrado hasta ahora
}
}
__END__
