sub
new {
my
(
$pkg
,
@args
)=
@_
;
my
$self
=
$pkg
->SUPER::new(
@args
);
my
(
$eNewick
,
$edgesR
,
$leavesR
,
$numleaves
,
$graph
,
$tree
,
$mudataR
)=
$self
->_rearrange([
qw(ENEWICK
EDGES
LEAVES
NUMLEAVES
GRAPH
TREE
MUDATA)
],
@args
);
bless
(
$self
,
$pkg
);
$self
->build_from_eNewick(
$eNewick
)
if
defined
$eNewick
;
$self
->build_from_edges(
@$edgesR
)
if
defined
$edgesR
;
$self
->build_from_graph(
$graph
)
if
defined
$graph
;
$self
->build_from_tree(
$tree
)
if
defined
$tree
;
if
((!
defined
$leavesR
) && (
defined
$numleaves
)) {
my
@leaves
=
map
{
"l$_"
} (1..
$numleaves
);
$leavesR
=\
@leaves
;
}
$self
->build_from_mudata(
$mudataR
,
$leavesR
)
if
((
defined
$mudataR
) && (
defined
$leavesR
));
return
$self
;
}
sub
build_from_edges {
my
(
$self
,
@edges
)=
@_
;
my
$graph
=Graph::Directed->new();
$graph
->add_edges(
@edges
);
$self
->{graph}=
$graph
;
$self
->recompute();
my
$labels
={};
foreach
my
$node
(
$self
->nodes()) {
$labels
->{
$node
}=
$node
;
}
$self
->{labels}=
$labels
;
}
sub
build_from_graph {
my
(
$self
,
$graph
)=
@_
;
my
$graphcp
=
$graph
->copy();
$self
->{graph}=
$graphcp
;
$self
->recompute();
my
$labels
={};
foreach
my
$node
(
$self
->nodes()) {
$labels
->{
$node
}=
$node
;
}
$self
->{labels}=
$labels
;
}
my
$_eN_index
;
sub
build_from_eNewick {
my
(
$self
,
$string
)=
@_
;
$_eN_index
=0;
my
$graph
=Graph::Directed->new();
my
$labels
={};
my
@blocks
=
split
(/; */,
$string
);
foreach
my
$block
(
@blocks
) {
my
(
$rt
,
$str
)=get_root_and_subtree(
$block
);
my
(
$rtlbl
,
$rttype
,
$rtid
,
$rtlng
)=get_label_type_id_length(
$rt
);
process_block(
$graph
,
$labels
,
$block
,
$rtid
);
$labels
->{
$rtid
}=
$rtlbl
.
''
;
}
$self
->{graph}=
$graph
;
$self
->{labels}=
$labels
;
$self
->recompute();
}
sub
process_block {
my
(
$graph
,
$labels
,
$block
,
$rtid
)=
@_
;
my
(
$rt
,
$str
)=get_root_and_subtree(
$block
);
my
@substrs
=my_split(
$str
);
foreach
my
$substr
(
@substrs
) {
my
(
$subrt
,
$subblock
)=get_root_and_subtree(
$substr
);
my
(
$subrtlbl
,
$subrttype
,
$subrtid
,
$subrtlng
)=
get_label_type_id_length(
$subrt
);
if
(!
$subrtlng
eq
''
) {
$graph
->add_weighted_edges(
$rtid
,
$subrtid
,
$subrtlng
);
}
else
{
$graph
->add_edges(
$rtid
,
$subrtid
);
}
if
(!
$subrttype
eq
''
) {
$graph
->set_edge_attribute(
$rtid
,
$subrtid
,
'type'
,
$subrttype
);
}
$subrtlbl
.=
''
;
if
((!
defined
$labels
->{
$subrtid
})||(
$labels
->{
$subrtid
} eq
''
)){
$labels
->{
$subrtid
}=
$subrtlbl
;
}
elsif
((
$labels
->{
$subrtid
} ne
$subrtlbl
)&&(
$subrtlbl
ne
''
)) {
die
(
"Different labels for the same node ("
.
$labels
->{
$subrtid
}.
" and $subrtlbl)"
);
}
if
(
$subblock
ne
""
) {
process_block(
$graph
,
$labels
,
$subblock
,
$subrtid
);
}
}
}
sub
get_root_and_subtree {
my
(
$block
)=
@_
;
my
(
$rt
,
$str
)=(
""
,
""
);
(
$rt
,
$str
)=
split
(/=/,
$block
);
if
(
$rt
eq
$block
) {
my
$pos
=
length
(
$rt
)-1;
while
((
substr
(
$rt
,
$pos
,1) ne
")"
) && (
$pos
>=0)) {
$pos
--;
}
$rt
=
substr
(
$block
,
$pos
+1,
length
(
$block
)-
$pos
);
$str
=
substr
(
$block
,0,
$pos
+1);
}
$rt
=trim(
$rt
);
$str
=trim(
$str
);
return
(
$rt
,
$str
);
}
sub
get_label_type_id_length {
my
(
$string
) =
@_
;
$string
.=
''
;
if
(
index
(
$string
,
'#'
)==-1) {
my
(
$label
,
$length
)=
split
(
':'
,
$string
);
$label
.=
''
;
my
$id
;
if
((!
defined
$label
) || (
$label
eq
''
)) {
$_eN_index
++;
$id
=
"T$_eN_index"
;
}
else
{
$id
=
$label
;
}
return
(
$label
,
''
,
$id
,
$length
);
}
else
{
my
(
$label
,
$string2
)=
split
(
'#'
,
$string
);
my
(
$typeid
,
$length
)=
split
(
':'
,
$string2
);
my
$type
=
$typeid
;
$type
=~ s/\d//g;
my
$id
=
$typeid
;
$id
=~ s/\D//g;
return
(
$label
,
$type
,
'#'
.
$id
,
$length
);
}
}
sub
trim
{
my
(
$string
) =
@_
;
$string
=~ s/^\s+//;
$string
=~ s/\s+$//;
return
$string
;
}
sub
my_split {
my
(
$string
) =
@_
;
my
$temp
=
""
;
my
@substrings
;
my
$level
=1;
for
my
$i
( 1 ..
length
(
$string
) ) {
my
$char
=
substr
(
$string
,
$i
,1);
if
(
$char
eq
"("
) {
$level
++;
}
if
(
$char
eq
")"
) {
if
(
$level
==1) {
push
@substrings
,
$temp
;
$temp
=
""
;
}
$level
--;
}
if
((
$char
eq
","
) && (
$level
==1)) {
push
@substrings
,
$temp
;
$temp
=
""
;
$char
=
""
;
}
$temp
=
$temp
.
$char
;
}
return
@substrings
;
}
sub
build_from_mudata {
my
(
$self
,
$mus
,
$leavesR
)=
@_
;
my
$graph
=Graph::Directed->new();
my
@nodes
=
keys
%{
$mus
};
my
@leaves
=@{
$leavesR
};
my
%seen
;
my
@internal
;
@seen
{
@leaves
} = ();
foreach
my
$node
(
@nodes
) {
push
(
@internal
,
$node
)
unless
exists
$seen
{
$node
};
}
@internal
=
sort
{
$mus
->{
$b
} <=>
$mus
->{
$a
} }
@internal
;
@nodes
=(
@internal
,
@leaves
);
my
$numnodes
=
@nodes
;
for
(
my
$i
=0;
$i
<
$numnodes
;
$i
++) {
my
$mu
=
$mus
->{
$nodes
[
$i
]};
my
$j
=
$i
+1;
while
(
$mu
->is_positive() &&
$j
<
$numnodes
) {
if
(
$mu
->geq_poset(
$mus
->{
$nodes
[
$j
]})) {
$graph
->add_edges((
$nodes
[
$i
],
$nodes
[
$j
]));
$mu
=
$mu
-
$mus
->{
$nodes
[
$j
]};
}
$j
++;
}
}
$self
->build_from_graph(
$graph
);
}
sub
build_from_tree {
my
(
$self
,
$tree
)=
@_
;
my
$str
;
my
$io
=IO::String->new(
$str
);
my
$treeio
=Bio::TreeIO->new(
-format
=>
'newick'
,
-fh
=>
$io
);
$treeio
->write_tree(
$tree
);
$self
->build_from_eNewick(
$str
);
}
sub
recompute {
my
(
$self
)=
@_
;
$self
->throw(
"Graph is not DAG:"
.
$self
->{graph})
unless
$self
->{graph}->is_dag();
my
@leaves
=
$self
->{graph}->successorless_vertices();
@leaves
=
sort
@leaves
;
my
$numleaves
=
@leaves
;
my
@roots
=
$self
->{graph}->predecessorless_vertices();
my
$numroots
=
@roots
;
my
@nodes
=
$self
->{graph}->vertices();
@nodes
=
sort
@nodes
;
my
$numnodes
=
@nodes
;
foreach
my
$node
(
@nodes
) {
if
(!
defined
$self
->{labels}->{
$node
}) {
$self
->{labels}->{
$node
}=
''
;
}
}
$self
->{leaves}=\
@leaves
;
$self
->{numleaves}=
$numleaves
;
$self
->{roots}=\
@roots
;
$self
->{numroots}=
$numroots
;
$self
->{nodes}=\
@nodes
;
$self
->{numnodes}=
$numnodes
;
$self
->{mudata}={};
$self
->{h}={};
$self
->compute_height();
$self
->compute_mu();
return
$self
;
}
sub
is_attackable {
my
(
$self
,
$u1
,
$v1
,
$u2
,
$v2
)=
@_
;
if
(
$self
->is_hybrid_node(
$v1
) ||
$self
->is_hybrid_node(
$v2
) ||
$self
->graph->is_reachable(
$v2
,
$u1
) ||
((
$u1
eq
$u2
)&&(
$v1
eq
$v2
)) ||
(!
scalar
grep
{(
$_
ne
$v2
) && (
$self
->is_tree_node(
$_
))}
$self
->graph->successors(
$u2
)))
{
return
0;
}
return
1;
}
sub
do_attack {
my
(
$self
,
$u1
,
$v1
,
$u2
,
$v2
,
$lbl
)=
@_
;
my
$graph
=
$self
->{graph};
$graph
->delete_edge(
$u1
,
$v1
);
$graph
->delete_edge(
$u2
,
$v2
);
$graph
->add_edge(
$u1
,
"T$lbl"
);
$graph
->add_edge(
"T$lbl"
,
$v1
);
$graph
->add_edge(
$u2
,
"#H$lbl"
);
$graph
->add_edge(
"#H$lbl"
,
$v2
);
$graph
->add_edge(
"T$lbl"
,
"#H$lbl"
);
$self
->build_from_graph(
$graph
);
}
sub
compute_mu {
my
(
$self
)=
@_
;
my
$graph
=
$self
->{graph};
my
$mudata
=
$self
->{mudata};
my
@leaves
=@{
$self
->{leaves}};
my
$numleaves
=
$self
->{numleaves};
for
(
my
$i
=0;
$i
<
$numleaves
;
$i
++) {
my
$vec
=Bio::PhyloNetwork::muVector->new(
$numleaves
);
$vec
->[
$i
]=1;
$mudata
->{
$leaves
[
$i
]}=
$vec
;
}
my
$h
=1;
while
(
my
@nodes
=
grep
{
$self
->{h}->{
$_
} ==
$h
} @{
$self
->{nodes}} )
{
foreach
my
$u
(
@nodes
) {
my
$vec
=Bio::PhyloNetwork::muVector->new(
$numleaves
);
foreach
my
$son
(
$graph
->successors(
$u
)) {
$vec
+=
$mudata
->{
$son
};
}
$mudata
->{
$u
}=
$vec
;
}
$h
++;
}
}
sub
compute_height {
my
(
$self
)=
@_
;
my
$graph
=
$self
->{graph};
my
@leaves
=@{
$self
->{leaves}};
foreach
my
$leaf
(
@leaves
) {
$self
->{h}->{
$leaf
}=0;
}
my
$h
=0;
while
(
my
@nodes
=
grep
{(
defined
$self
->{h}->{
$_
})&&(
$self
->{h}->{
$_
} ==
$h
)}
@{
$self
->{nodes}} )
{
foreach
my
$node
(
@nodes
) {
foreach
my
$parent
(
$graph
->predecessors(
$node
)) {
$self
->{h}->{
$parent
}=
$h
+1;
}
}
$h
++;
}
}
sub
is_leaf {
my
(
$self
,
$node
)=
@_
;
if
(
$self
->{graph}->out_degree(
$node
) == 0) {
return
1;}
return
0;
}
sub
is_root {
my
(
$self
,
$node
)=
@_
;
if
(
$self
->{graph}->in_degree(
$node
) == 0) {
return
1;}
return
0;
}
sub
is_tree_node {
my
(
$self
,
$node
)=
@_
;
if
(
$self
->{graph}->in_degree(
$node
) <= 1) {
return
1;}
return
0;
}
sub
is_hybrid_node {
my
(
$self
,
$node
)=
@_
;
if
(
$self
->{graph}->in_degree(
$node
) > 1) {
return
1;}
return
0;
}
sub
has_tree_child {
my
$g
=
shift
(
@_
);
my
$node
=
shift
(
@_
);
my
@Sons
=
$g
->successors(
$node
);
foreach
my
$son
(
@Sons
) {
if
(
$g
->in_degree(
$son
)==1) {
return
1;
}
}
return
0;
}
sub
is_tree_child {
my
(
$self
)=
@_
;
if
(
defined
$self
->{is_tree_child}) {
return
$self
->{is_tree_child};
}
$self
->{is_tree_child}=0;
my
$graph
=
$self
->{graph};
foreach
my
$node
(@{
$self
->{nodes}}) {
return
0
unless
(
$graph
->out_degree(
$node
)==0 ||
has_tree_child(
$graph
,
$node
));
}
$self
->{is_tree_child}=1;
return
1;
}
sub
nodes {
my
(
$self
)=
@_
;
return
@{
$self
->{nodes}};
}
sub
leaves {
my
(
$self
)=
@_
;
return
@{
$self
->{leaves}};
}
sub
roots {
my
(
$self
)=
@_
;
return
@{
$self
->{roots}};
}
sub
internal_nodes {
my
(
$self
)=
@_
;
return
grep
{!
$self
->is_leaf(
$_
)}
$self
->nodes();
}
sub
tree_nodes {
my
(
$self
)=
@_
;
return
grep
{
$self
->is_tree_node(
$_
)}
$self
->nodes();
}
sub
hybrid_nodes {
my
(
$self
)=
@_
;
return
grep
{
$self
->is_hybrid_node(
$_
)}
$self
->nodes();
}
sub
graph {
my
(
$self
)=
@_
;
return
$self
->{graph};
}
sub
edges {
my
(
$self
)=
@_
;
return
$self
->{graph}->edges();
}
sub
tree_edges {
my
(
$self
)=
@_
;
return
grep
{
$self
->is_tree_node(
$_
->[1])}
$self
->edges();
}
sub
hybrid_edges {
my
(
$self
)=
@_
;
return
grep
{
$self
->is_hybrid_node(
$_
->[1])}
$self
->edges();
}
sub
explode {
my
(
$self
)=
@_
;
my
@trees
;
$self
->explode_rec(\
@trees
);
return
@trees
;
}
sub
explode_rec {
my
(
$self
,
$trees
)=
@_
;
my
@h
=
$self
->hybrid_nodes;
if
(
scalar
@h
) {
my
$v
=
shift
@h
;
for
my
$u
(
$self
->{graph}->predecessors(
$v
)) {
$self
->{graph}->delete_edge(
$u
,
$v
);
$self
->explode_rec(
$trees
);
$self
->{graph}->add_edge(
$u
,
$v
);
}
}
else
{
my
$io
= IO::String->new(
$self
->eNewick);
my
$treeio
= Bio::TreeIO->new(
-format
=>
'newick'
,
-fh
=>
$io
);
my
$tree
=
$treeio
->next_tree;
$tree
->contract_linear_paths;
push
@{
$trees
},
$tree
;
}
}
sub
mudata {
my
(
$self
)=
@_
;
return
%{
$self
->{mudata}};
}
sub
mudata_node {
my
(
$self
,
$u
)=
@_
;
return
$self
->{mudata}{
$u
};
}
sub
heights {
my
(
$self
)=
@_
;
return
%{
$self
->{h}};
}
sub
height_node {
my
(
$self
,
$u
)=
@_
;
return
$self
->{h}{
$u
};
}
sub
mu_distance {
my
(
$net1
,
$net2
)=
@_
;
my
@nodes1
=
$net1
->nodes;
my
@nodes2
=
$net2
->nodes;
my
$comp
= Array::Compare->new;
$net1
->throw(
"Cannot compare phylogenetic networks on different set of leaves"
)
unless
$comp
->compare(
$net1
->{leaves},
$net2
->{leaves});
$net1
->
warn
(
"Not a tree-child phylogenetic network"
)
unless
$net1
->is_tree_child();
$net2
->
warn
(
"Not a tree-child phylogenetic network"
)
unless
$net2
->is_tree_child();
my
@leaves
=@{
$net1
->{leaves}};
my
%matched1
;
my
%matched2
;
OUTER:
foreach
my
$node1
(
@nodes1
) {
foreach
my
$node2
(
@nodes2
) {
if
(
(!
exists
$matched1
{
$node1
}) && (!
exists
$matched2
{
$node2
}) &&
(
$net1
->{mudata}{
$node1
} ==
$net2
->{mudata}{
$node2
})
) {
$matched1
{
$node1
}=
$node2
;
$matched2
{
$node2
}=
$node1
;
next
OUTER;
}
}
}
return
(
scalar
@nodes1
)+(
scalar
@nodes2
)-2*(
scalar
keys
%matched1
);
}
sub
mu_distance_generalized {
my
(
$net1
,
$net2
)=
@_
;
my
(
$netr1
,
$netr2
)=
$net1
->topological_restriction(
$net2
);
return
$netr1
->mu_distance(
$netr2
);
}
sub
mudata_string_node {
my
(
$self
,
$u
)=
@_
;
return
$self
->{mudata}->{
$u
}->display();
}
sub
mudata_string {
my
(
$self
)=
@_
;
return
$self
->{mudata_string}
if
defined
$self
->{mudata_string};
my
@internal
=
$self
->internal_nodes;
my
$mus
=
$self
->{mudata};
@internal
=
sort
{
$mus
->{
$b
} <=>
$mus
->{
$a
} }
@internal
;
my
$str
=
""
;
foreach
my
$node
(
@internal
) {
$str
=
$str
.
$self
->mudata_string_node(
$node
);
}
$self
->{mudata_string}=
$str
;
return
$str
;
}
sub
is_mu_isomorphic {
my
(
$net1
,
$net2
)=
@_
;
return
(
$net1
->mudata_string() eq
$net2
->mudata_string());
}
sub
compute_tripartition_node {
my
(
$self
,
$u
)=
@_
;
$self
->
warn
(
"Cannot compute tripartitions on unrooted networks. Will assume one at random"
)
unless
(
$self
->{numroots} == 1);
my
$root
=
$self
->{roots}->[0];
my
$graph
=
$self
->{graph};
my
$graphPruned
=
$graph
->copy();
$graphPruned
->delete_vertex(
$u
);
my
$tripartition
=
""
;
foreach
my
$leaf
(@{
$self
->{leaves}}) {
my
$type
;
if
(
$graph
->is_reachable(
$u
,
$leaf
)) {
if
(
$graphPruned
->is_reachable(
$root
,
$leaf
)) {
$type
=
"B"
;}
else
{
$type
=
"A"
;}
}
else
{
$type
=
"C"
;}
$tripartition
.=
$type
;
}
$self
->{tripartitions}->{
$u
}=
$tripartition
;
}
sub
compute_tripartitions {
my
(
$self
)=
@_
;
foreach
my
$node
(@{
$self
->{nodes}}) {
$self
->compute_tripartition_node(
$node
);
}
}
sub
tripartitions {
my
(
$self
)=
@_
;
$self
->compute_tripartitions()
unless
defined
$self
->{tripartitions};
return
%{
$self
->{tripartitions}};
}
sub
tripartition_error {
my
(
$net1
,
$net2
)=
@_
;
my
$comp
= Array::Compare->new;
$net1
->throw(
"Cannot compare phylogenetic networks on different set of leaves"
)
unless
$comp
->compare(
$net1
->{leaves},
$net2
->{leaves});
$net1
->
warn
(
"Not a tree-child phylogenetic network"
)
unless
$net1
->is_tree_child();
$net2
->
warn
(
"Not a tree-child phylogenetic network"
)
unless
$net2
->is_tree_child();
$net1
->
warn
(
"Not a time-consistent network"
)
unless
$net1
->is_time_consistent();
$net2
->
warn
(
"Not a time-consistent network"
)
unless
$net2
->is_time_consistent();
$net1
->compute_tripartitions()
unless
defined
$net1
->{tripartitions};
$net2
->compute_tripartitions()
unless
defined
$net2
->{tripartitions};
my
@edges1
=
$net1
->{graph}->edges();
my
@edges2
=
$net2
->{graph}->edges();
my
(
$FN
,
$FP
)=(0,0);
foreach
my
$edge1
(
@edges1
) {
my
$matched
=0;
foreach
my
$edge2
(
@edges2
) {
if
(
$net1
->{tripartitions}->{
$edge1
->[1]} eq
$net2
->{tripartitions}->{
$edge2
->[1]}) {
$matched
=1;
last
;
}
}
if
(!
$matched
) {
$FN
++;}
}
foreach
my
$edge2
(
@edges2
) {
my
$matched
=0;
foreach
my
$edge1
(
@edges1
) {
if
(
$net1
->{tripartitions}->{
$edge1
->[1]} eq
$net2
->{tripartitions}->{
$edge2
->[1]}) {
$matched
=1;
last
;
}
}
if
(!
$matched
) {
$FP
++;}
}
return
(
$FN
/(
scalar
@edges1
)+
$FP
/(
scalar
@edges2
))/2;
}
sub
is_time_consistent {
my
(
$self
)=
@_
;
$self
->compute_temporal_representation()
unless
exists
$self
->{has_temporal_representation};
return
$self
->{has_temporal_representation};
}
sub
temporal_representation {
my
(
$self
)=
@_
;
if
(
$self
->is_time_consistent) {
return
%{
$self
->{temporal_representation}};
}
return
0;
}
sub
compute_temporal_representation {
my
(
$self
)=
@_
;
my
$quotient
=Graph::Directed->new();
my
$classes
=find_classes(
$self
);
my
%repr
;
map
{
$repr
{
$_
}=
$classes
->{
$_
}[0]}
$self
->nodes();
foreach
my
$e
(
$self
->tree_edges()) {
$quotient
->add_edge(
$repr
{
$e
->[0]},
$repr
{
$e
->[1]});
}
my
%temp
;
my
$depth
=0;
while
(
$quotient
->vertices()) {
if
(
my
@svs
=
$quotient
->predecessorless_vertices()) {
foreach
my
$sv
(
@svs
) {
$temp
{
$sv
}=
$depth
;
}
$quotient
->delete_vertices(
@svs
);
}
else
{
return
0;
}
$depth
++;
}
foreach
my
$node
(@{
$self
->{nodes}}) {
$temp
{
$node
}=
$temp
{
$repr
{
$node
}}
}
$self
->{temporal_representation}=\
%temp
;
$self
->{has_temporal_representation}=1;
}
sub
find_classes {
my
(
$self
)=
@_
;
my
$classes
={};
map
{
$classes
->{
$_
}=[
$_
]}
$self
->nodes();
foreach
my
$e
(
$self
->hybrid_edges()) {
$classes
=join_classes(
$classes
,
$e
->[0],
$e
->[1]);
}
return
$classes
;
}
sub
join_classes {
my
(
$classes
,
$u
,
$v
)=
@_
;
my
@clu
=@{
$classes
->{
$u
}};
my
@clv
=@{
$classes
->{
$v
}};
my
@cljoin
=(
@clu
,
@clv
);
map
{
$classes
->{
$_
}=\
@cljoin
}
@cljoin
;
return
$classes
;
}
sub
contract_elementary {
my
(
$self
)=
@_
;
my
$contracted
=
$self
->graph->copy();
my
@nodes
=
$self
->nodes();
my
$mus
=
$self
->{mudata};
my
$hs
=
$self
->{h};
my
%blocks
;
foreach
my
$u
(
@nodes
) {
$blocks
{
$u
}=[
$u
];
}
my
@elementary
=
grep
{
$contracted
->out_degree(
$_
) == 1}
$self
->tree_nodes();
@elementary
=
sort
{
$mus
->{
$b
} <=>
$mus
->{
$a
} ||
$hs
->{
$b
} <=>
$hs
->{
$a
}}
@elementary
;
foreach
my
$elem
(
@elementary
) {
my
@children
=
$contracted
->successors(
$elem
);
my
$child
=
$children
[0];
if
(
$contracted
->in_degree(
$elem
) == 1) {
my
@parents
=
$contracted
->predecessors(
$elem
);
my
$parent
=
$parents
[0];
$contracted
->add_edge(
$parent
,
$child
);
}
$contracted
->delete_vertex(
$elem
);
my
@blch
=@{
$blocks
{
$child
}};
my
@blem
=@{
$blocks
{
$elem
}};
$blocks
{
$child
}=[
@blem
,
@blch
];
delete
$blocks
{
$elem
};
}
my
$contr
=Bio::PhyloNetwork->new(
-graph
=>
$contracted
);
return
$contr
,\
%blocks
;
}
sub
optimal_alignment {
my
(
$net1
,
$net2
,
%params
)=
@_
;
my
(
$net1cont
,
$blocks1
)=contract_elementary(
$net1
);
my
(
$net2cont
,
$blocks2
)=contract_elementary(
$net2
);
my
(
$wc
,
$alignc
,
$weightc
)=
optimal_alignment_noelementary(
$net1cont
,
$net2cont
,
%params
);
my
%alignment
=();
my
$totalweigth
=0;
my
%weigths
=();
foreach
my
$u1
(
keys
%$alignc
) {
my
$u2
=
$alignc
->{
$u1
};
my
@block1
=@{
$blocks1
->{
$u1
}};
my
@block2
=@{
$blocks2
->{
$u2
}};
while
(
@block1
&&
@block2
) {
my
$u1dc
=
pop
@block1
;
my
$u2dc
=
pop
@block2
;
$alignment
{
$u1dc
}=
$u2dc
;
$weigths
{
$u1dc
}=
$weightc
->{
$u1
};
$totalweigth
+=
$weigths
{
$u1dc
};
}
}
return
$totalweigth
,\
%alignment
,\
%weigths
;
}
sub
optimal_alignment_noelementary {
my
(
$net1
,
$net2
,
%params
)=
@_
;
my
$comp
= Array::Compare->new;
$net1
->throw(
"Cannot align phylogenetic networks on different set of leaves"
)
unless
$comp
->compare(
$net1
->{leaves},
$net2
->{leaves});
my
$distance
;
if
((
defined
$params
{-metric})and (
$params
{-metric} eq
'Hamming'
)) {
$distance
=
'Hamming'
;
}
else
{
$distance
=
'Manhattan'
;
}
my
$numleaves
=
$net1
->{numleaves};
my
@nodes1
=
$net1
->internal_nodes();
my
@nodes2
=
$net2
->internal_nodes();
my
$numnodes1
=
@nodes1
;
my
$numnodes2
=
@nodes2
;
my
@matrix
=();
for
(
my
$i
=0;
$i
<
$numnodes1
;
$i
++) {
my
@row
=();
for
(
my
$j
=0;
$j
<
$numnodes2
;
$j
++) {
push
@row
,weight(
$net1
,
$nodes1
[
$i
],
$net2
,
$nodes2
[
$j
],
$distance
);
}
push
@matrix
,\
@row
;
}
my
@alignment
=();
Algorithm::Munkres::assign(\
@matrix
,\
@alignment
);
my
%alignmenthash
;
my
%weighthash
;
my
$totalw
=0;
foreach
my
$leaf
(@{
$net1
->{leaves}}) {
$alignmenthash
{
$leaf
}=
$leaf
;
$weighthash
{
$leaf
}=0;
}
for
(
my
$i
=0;
$i
<
$numnodes1
;
$i
++) {
if
(
defined
$nodes2
[
$alignment
[
$i
]]) {
$alignmenthash
{
$nodes1
[
$i
]}=
$nodes2
[
$alignment
[
$i
]];
$weighthash
{
$nodes1
[
$i
]}=
$matrix
[
$i
][
$alignment
[
$i
]];
$totalw
+=
$matrix
[
$i
][
$alignment
[
$i
]];
}
}
return
$totalw
,\
%alignmenthash
,\
%weighthash
;
}
sub
optimal_alignment_generalized {
my
(
$net1
,
$net2
,
%params
)=
@_
;
my
(
$netr1
,
$netr2
)=
$net1
->topological_restriction(
$net2
);
return
$netr1
->optimal_alignment(
$netr2
,
%params
);
}
sub
weight {
my
(
$net1
,
$v1
,
$net2
,
$v2
,
$distance
)=
@_
;
my
$w
;
if
(!
defined
$distance
) {
$distance
=
'Manhattan'
;
}
if
(
$distance
eq
'Hamming'
) {
$w
=
$net1
->{mudata}->{
$v1
}->hamming(
$net2
->{mudata}->{
$v2
});
}
else
{
$w
=
$net1
->{mudata}->{
$v1
}->manhattan(
$net2
->{mudata}->{
$v2
});
}
if
((
$net1
->is_tree_node(
$v1
) &&
$net2
->is_hybrid_node(
$v2
)) ||
(
$net2
->is_tree_node(
$v2
) &&
$net1
->is_hybrid_node(
$v1
))
)
{
$w
+=1/(2
*$net1
->{numleaves});
}
return
$w
;
}
sub
topological_restriction {
my
(
$net1
,
$net2
)=
@_
;
my
@leaves1
=
$net1
->leaves();
my
@leaves2
=
$net2
->leaves();
my
$numleaves1
=
scalar
@leaves1
;
my
$numleaves2
=
scalar
@leaves2
;
my
%position1
;
for
(
my
$i
=0;
$i
<
$numleaves1
;
$i
++) {
$position1
{
$leaves1
[
$i
]}=
$i
;
}
my
%position2
;
my
@commonleaves
=();
for
(
my
$j
=0;
$j
<
$numleaves2
;
$j
++) {
if
(
defined
$position1
{
$leaves2
[
$j
]}) {
push
@commonleaves
,
$leaves2
[
$j
];
$position2
{
$leaves2
[
$j
]}=
$j
;
}
}
my
$graphred1
=
$net1
->{graph}->copy();
my
$graphred2
=
$net2
->{graph}->copy();
OUTER1:
foreach
my
$u
(
$graphred1
->vertices()) {
my
$mu
=
$net1
->mudata_node(
$u
);
foreach
my
$leaf
(
@commonleaves
) {
if
(
$mu
->[
$position1
{
$leaf
}]>0) {
next
OUTER1;
}
}
$graphred1
->delete_vertex(
$u
);
}
OUTER2:
foreach
my
$u
(
$graphred2
->vertices()) {
my
$mu
=
$net2
->mudata_node(
$u
);
foreach
my
$leaf
(
@commonleaves
) {
if
(
$mu
->[
$position2
{
$leaf
}]>0) {
next
OUTER2;
}
}
$graphred2
->delete_vertex(
$u
);
}
my
$netr1
=Bio::PhyloNetwork->new(
-graph
=>
$graphred1
);
my
$netr2
=Bio::PhyloNetwork->new(
-graph
=>
$graphred2
);
return
(
$netr1
,
$netr2
);
}
sub
eNewick {
my
(
$self
)=
@_
;
my
$str
=
""
;
my
$seen
={};
foreach
my
$root
(
$self
->roots()) {
$str
=
$str
.
$self
->eNewick_aux(
$root
,
$seen
,
undef
).
"; "
;
}
return
$str
;
}
sub
eNewick_aux {
my
(
$self
,
$node
,
$seen
,
$parent
)=
@_
;
my
$str
=
''
;
if
(
$self
->is_leaf(
$node
) ||
(
defined
$seen
->{
$node
}) )
{
$str
=make_label(
$self
,
$parent
,
$node
);
}
else
{
$seen
->{
$node
}=1;
my
@sons
=
$self
->{graph}->successors(
$node
);
$str
=
"("
;
foreach
my
$son
(
@sons
) {
$str
=
$str
.
$self
->eNewick_aux(
$son
,
$seen
,
$node
).
","
;
}
chop
(
$str
);
$str
.=
")"
.make_label(
$self
,
$parent
,
$node
);
}
return
$str
;
}
sub
make_label {
my
(
$self
,
$parent
,
$node
)=
@_
;
my
$str
=
''
;
if
(
$self
->is_hybrid_node(
$node
)) {
my
$lbl
=
$self
->{labels}->{
$node
};
if
(
$lbl
=~ /
$lbl
=
''
;
}
$str
.=
$lbl
;
$str
.=
'#'
;
if
((
defined
$parent
) &&
(
$self
->graph->has_edge_attribute(
$parent
,
$node
,
'type'
))) {
$str
.=
$self
->graph->get_edge_attribute(
$parent
,
$node
,
'type'
);
}
$str
.=
substr
$node
,1;
}
else
{
$str
.=
$self
->{labels}->{
$node
};
}
if
((
defined
$parent
) &&
(
$self
->graph->has_edge_weight(
$parent
,
$node
))) {
$str
.=
":"
.
$self
->graph->get_edge_weight(
$parent
,
$node
);
}
return
$str
;
}
sub
eNewick_full {
my
(
$self
)=
@_
;
my
$str
=
""
;
my
$seen
={};
foreach
my
$root
(
$self
->roots()) {
$str
=
$str
.
$self
->eNewick_full_aux(
$root
,
$seen
,
undef
).
"; "
;
}
return
$str
;
}
sub
eNewick_full_aux {
my
(
$self
,
$node
,
$seen
,
$parent
)=
@_
;
my
$str
=
''
;
if
(
$self
->is_leaf(
$node
) ||
(
defined
$seen
->{
$node
}) )
{
$str
=make_label_full(
$self
,
$parent
,
$node
);
}
else
{
$seen
->{
$node
}=1;
my
@sons
=
$self
->{graph}->successors(
$node
);
$str
=
"("
;
foreach
my
$son
(
@sons
) {
$str
=
$str
.
$self
->eNewick_full_aux(
$son
,
$seen
,
$node
).
","
;
}
chop
(
$str
);
$str
.=
")"
.make_label_full(
$self
,
$parent
,
$node
);
}
return
$str
;
}
sub
make_label_full {
my
(
$self
,
$parent
,
$node
)=
@_
;
my
$str
=
''
;
if
(
$self
->is_hybrid_node(
$node
)) {
my
$lbl
=
$self
->{labels}->{
$node
};
if
(
$lbl
=~ /
$lbl
=
''
;
}
$str
.=
$lbl
;
$str
.=
'#'
;
if
((
defined
$parent
) &&
(
$self
->graph->has_edge_attribute(
$parent
,
$node
,
'type'
))) {
$str
.=
$self
->graph->get_edge_attribute(
$parent
,
$node
,
'type'
);
}
$str
.=
substr
$node
,1;
}
else
{
if
((
defined
$self
->{labels}->{
$node
})&&(
$self
->{labels}->{
$node
} ne
''
)) {
$str
.=
$self
->{labels}->{
$node
};
}
else
{
$str
.=
$node
;
}
}
if
((
defined
$parent
) &&
(
$self
->graph->has_edge_weight(
$parent
,
$node
))) {
$str
.=
":"
.
$self
->graph->get_edge_weight(
$parent
,
$node
);
}
return
$str
;
}
sub
display {
my
(
$self
)=
@_
;
my
$str
=
""
;
my
$graph
=
$self
->{graph};
my
@leaves
=
$self
->leaves();
my
@nodes
=@{
$self
->{nodes}};
$str
.=
"Leaves:\t@leaves\n"
;
$str
.=
"Nodes:\t@nodes\n"
;
$str
.=
"Graph:\t$graph\n"
;
$str
.=
"eNewick:\t"
.
$self
->eNewick().
"\n"
;
$str
.=
"Full eNewick:\t"
.
$self
->eNewick_full().
"\n"
;
$str
.=
"Mu-data and heights:\n"
;
foreach
my
$node
(
@nodes
) {
$str
.=
"v=$node: "
;
if
(
exists
$self
->{labels}->{
$node
}) {
$str
.=
"\tlabel="
.
$self
->{labels}->{
$node
}.
","
;
}
else
{
$str
.=
"\tlabel=(none),"
;
}
$str
.=
"\th="
.
$self
->{h}->{
$node
}.
", \tmu="
.
$self
->{mudata}->{
$node
}.
"\n"
;
}
if
(
exists
$self
->{has_temporal_representation}) {
$str
.=
"Temporal representation:\n"
;
if
(
$self
->{has_temporal_representation}) {
foreach
my
$node
(
@nodes
) {
$str
.=
"v=$node; "
;
$str
.=
"\tt="
.
$self
->{temporal_representation}->{
$node
}.
"\n"
;
}
}
else
{
$str
.=
"Does not exist.\n"
;
}
}
if
(
exists
$self
->{tripartitions}) {
$str
.=
"Tripartitions:\n"
;
foreach
my
$node
(
@nodes
) {
$str
.=
"v=$node; "
;
$str
.=
"\ttheta="
.
$self
->{tripartitions}->{
$node
}.
"\n"
;
}
}
return
$str
;
}
1;