<html><head><meta name="color-scheme" content="light dark"></head><body><pre style="word-wrap: break-word; white-space: pre-wrap;">#!/usr/bin/perl -I .
# Author: Debashis Sahoo &lt;dsahoo@ucsd.edu&gt;

if (scalar(@ARGV) &lt;= 0) {
    print "perl analyze.pl &lt;cmd&gt; &lt;args&gt;\n";
    exit(1);
}

use U;
use Network;
use Graph;
use Hegemon;
use LWP::UserAgent;
use Data::Dumper;
use JSON;

my $cmd = shift(@ARGV);
if ($cmd eq "toidx") {
  &amp;U::convertidx(@ARGV);
}
if ($cmd eq "thr") {
  &amp;U::printThr(@ARGV);
}
if ($cmd eq "Info") {
  &amp;convertInfo(@ARGV);
}
if ($cmd eq "VInfo") {
  &amp;convertVInfo(@ARGV);
}
if ($cmd eq "bv") {
  &amp;convertBv(@ARGV);
}
if ($cmd eq "data-download") {
  &amp;downloadData(@ARGV);
}
if ($cmd eq "ibd") {
  &amp;analyzeIBD(@ARGV);
}

sub convertInfo {
  my ($dataset, $thr, @rest) = @_;
  my $pre = $dataset;
  $thr = "$pre\-thr.txt" if (!defined $thr);
  my $h = Hegemon-&gt;new(expr =&gt; "$pre\-expr.txt",
      idx =&gt; "$pre\-idx.txt", thr =&gt; $thr,
      survival =&gt; undef);
  $h-&gt;printInfo();
}

sub convertVInfo {
  my ($dataset, $thr, @rest) = @_;
  my $pre = $dataset;
  $thr = "$pre\-thr.txt" if (!defined $thr);
  my $h = Hegemon-&gt;new(expr =&gt; "$pre\-expr.txt",
      idx =&gt; "$pre\-idx.txt", thr =&gt; $thr,
      survival =&gt; undef);
  $h-&gt;printVInfo();
}

sub convertBv {
  my ($dataset, @rest) = @_;
  my $pre = $dataset;
  my $h = Hegemon-&gt;new(expr =&gt; "$pre\-expr.txt",
      idx =&gt; "$pre\-idx.txt", thr =&gt; "$pre\-thr.txt",
      survival =&gt; undef);
  $h-&gt;printBvFile();
}

sub gamma {
  my ($u, $edges, $hash) = @_;
  if (defined $hash &amp;&amp; defined $hash-&gt;{$u}) {
    return $hash-&gt;{$u};
  }
  my $res = [$u, keys %{$edges-&gt;{$u}}];
  if (defined $hash) {
    $hash-&gt;{$u} = $res;
  }
  return $res;
}
sub rho {
  my ($u, $v, $edges, $hash, $shash) = @_;
  if (defined $shash &amp;&amp; defined $shash-&gt;{$u}-&gt;{$v}) {
    return $shash-&gt;{$u}-&gt;{$v};
  }
  my $gu = &amp;gamma($u, $edges, $hash);
  my $gv = &amp;gamma($v, $edges, $hash);
  my $g_union = &amp;U::union($gu, $gv);
  my $g_int = &amp;U::intersection($gu, $gv);
  my $res = scalar(@$g_int) / scalar(@$g_union);
  if (defined $shash) {
    $shash-&gt;{$u}-&gt;{$v} = $res;
  }
  return $res;
}
#https://arxiv.org/pdf/1606.00950.pdf
#https://link.springer.com/chapter/10.1007/978-3-319-91452-7_13
sub DCT {
  my ($nodes, $edges) = @_;
  my $T = [];
  my $S = [];
  my $V = [keys %{$nodes}];
  my $VB = {};
  my $checked = {};
  my $density = {};
  my $connect = {};
  my $count = 0;
  my $num = scalar(@$V);
  if (scalar(@$V) &gt; 0) {
    foreach my $u (@$V) {
      if (scalar(keys %{$edges-&gt;{$u}}) == 0) {
        $checked-&gt;{$u} = 1;
        push @$S, $u;
        $num --;
      }
      else {
        $VB-&gt;{$u} = 1;
      }
    }
    foreach my $u (@$V) {
      if (!defined $checked-&gt;{$u} &amp;&amp; defined $VB-&gt;{$u}) {
        $checked-&gt;{$u} = 1;
        push @$T, $u;
        $count ++;
        last;
      }
    }
  }
  my $hash = {};
  while ($count &lt; $num) {
    my $maxv = -1;
    my ($p, $q);
    foreach my $u (@$T) {
      my $gu = &amp;gamma($u, $edges, $hash);
      foreach my $v (@$gu) {
        if (!defined $checked-&gt;{$v} &amp;&amp; defined $VB-&gt;{$v}) {
          my $r = &amp;rho($u, $v, $edges, $hash);
          if ($r &gt; $maxv) {
            $maxv = $r;
            $p = $v;
            $q = $u;
          }
        }
      }
    }
    if (defined $p) {
      $checked-&gt;{$p} = 1;
      $connect-&gt;{$p} = $q;
      $density-&gt;{$p} = $maxv;
      push @$T, $p;
      $count++;
    }
    else {
      foreach my $u (@$V) {
        if (!defined $checked-&gt;{$u} &amp;&amp; defined $VB-&gt;{$u}) {
          $checked-&gt;{$u} = 1;
          push @$T, $u;
          $count++;
          $p = $u;
          my $gu = &amp;gamma($u, $edges);
          $maxv = scalar(@$gu);
          last;
        }
      }
    }
    print STDERR "$count/$num $p $q $maxv\n";
  }
  return [$T, $density, $connect];
}

sub DCT1 {
  my ($nodes, $edges) = @_;
  print STDERR scalar(keys %{$nodes}), "\n";
  my $V = [keys %{$nodes}];
  my $VB = {};
  foreach my $u (@$V) {
    if (scalar(keys %{$edges-&gt;{$u}}) &gt; 0) {
      $VB-&gt;{$u} = 1;
    }
  }
  my ($nodes_2, $edges_2) = &amp;Graph::subGraph1($nodes, $edges, $VB);
  print STDERR scalar(keys %{$nodes_2}), "\n";
  return &amp;DCT2($nodes_2, $edges_2);
}

sub DCT2 {
  my ($nodes, $edges) = @_;
  my $T = [];
  my $V = [keys %{$nodes}];
  my $checked = {};
  my $density = {};
  my $connect = {};
  my $count = 0;
  my $num = scalar(@$V);
  if (scalar(@$V) &gt; 0) {
    foreach my $u (@$V) {
      if (!defined $checked-&gt;{$u}) {
        $checked-&gt;{$u} = 1;
        push @$T, $u;
        $count ++;
        last;
      }
    }
  }
  my $hash = {};
  my $shash = {};
  while ($count &lt; $num) {
    my $maxv = -1;
    my ($p, $q);
    foreach my $u ($T-&gt;[$#{$T}]) {
      my $gu = &amp;gamma($u, $edges);
      foreach my $v (@$gu) {
        if (!defined $checked-&gt;{$v}) {
          my $r = &amp;rho($u, $v, $edges);
          if ($r &gt; $maxv) {
            $maxv = $r;
            $p = $v;
            $q = $u;
          }
        }
      }
    }
    if (defined $p) {
      $checked-&gt;{$p} = 1;
      $connect-&gt;{$p} = $q;
      $density-&gt;{$p} = $maxv;
      push @$T, $p;
      $count++;
    }
    else {
      foreach my $u (@$V) {
        if (!defined $checked-&gt;{$u}) {
          $checked-&gt;{$u} = 1;
          push @$T, $u;
          $count++;
          $p = $u;
          my $gu = &amp;gamma($u, $edges, $hash);
          $maxv = scalar(@$gu);
          last;
        }
      }
    }
    print STDERR "$count/$num $p $q $maxv\n";
  }
  return [$T, $density, $connect];
}

sub getCodes {
  my ($n, $l, $balanced, $balancedhash) = @_;
  my $num = $n-&gt;getNum();
  my $relations = {};
  foreach my $u (@$l) {
    my $i = $balancedhash-&gt;{$u};
    for (my $j = 0; $j &lt; $num/2; $j++) {
      my $code = $n-&gt;readCode_1_1($i, $j);
      if ($code &gt; 0) {
        if (!defined $relations-&gt;{$code}) {
          $relations-&gt;{$code} = {};
        }
        $relations-&gt;{$code}-&gt;{$balanced-&gt;[$j]} = 1;
      }
    }
  }
  return $relations;
}

sub convertCodes {
  my $vcode = shift;
  my $b0 = $vcode-&gt;bit_test(1);
  my $b1 = $vcode-&gt;bit_test(0);
  my $b2 = $vcode-&gt;bit_test(3);
  my $b3 = $vcode-&gt;bit_test(2);
  my $total = $b0 + $b1 + $b2 + $b3;
  if ($total == 1) {
    if ($b0 == 1) { return 1; }
    if ($b1 == 1) { return 2; }
    if ($b2 == 1) { return 3; }
    if ($b3 == 1) { return 4; }
  }
  if ($total == 2) {
    if ($b1 == 1 &amp;&amp; $b2 == 1) { return 5; }
    if ($b0 == 1 &amp;&amp; $b3 == 1) { return 6; }
  }
  return 0;
}

sub getNCodes {
  my ($n, $net, $l, $balanced, $balancedhash) = @_;
  my $num = $n-&gt;getNum();
  my $relations = {};
  my $nb = $n-&gt;{'numBits'};
  my $vcode = Bit::Vector-&gt;new($nb*2*2);
  foreach my $u (@$l) {
    my $i = $balancedhash-&gt;{$u};
    for (my $j = 0; $j &lt; $num/2; $j++) {
      $n-&gt;readNCode($net, $i, $j, $vcode);
      my $code = &amp;convertCodes($vcode);
      if ($code &gt; 0) {
        if (!defined $relations-&gt;{$code}) {
          $relations-&gt;{$code} = {};
        }
        $relations-&gt;{$code}-&gt;{$balanced-&gt;[$j]}++;
      }
    }
  }
  return $relations;
}

sub getNCodesArijs {
  my ($arijs, $n, $net, $l, $balanced, $balancedhash) = @_;
  my $num = $n-&gt;getNum();
  my $relations = {};
  my $nb = $n-&gt;{'numBits'};
  my $vcode = Bit::Vector-&gt;new($nb*2*2);
  foreach my $u (@$l) {
    my $i = $balancedhash-&gt;{$u};
    for (my $j = 0; $j &lt; $num/2; $j++) {
      $n-&gt;readNCode($net, $i, $j, $vcode);
      my $code = &amp;convertCodes($vcode);
      next if ($code &lt;= 0);
      my $res = &amp;getArijsCode($arijs, $u, $balanced-&gt;[$j]);
      if (&amp;isCodeCompatible($code, $res)) {
        if (!defined $relations-&gt;{$code}) {
          $relations-&gt;{$code} = {};
        }
        $relations-&gt;{$code}-&gt;{$balanced-&gt;[$j]}++;
      }
    }
  }
  return $relations;
}

sub getNCodesPair {
  my ($n, $net, $u, $v, $balanced, $balancedhash) = @_;
  my $num = $n-&gt;getNum();
  my $relations = {};
  my $nb = $n-&gt;{'numBits'};
  my $vcode = Bit::Vector-&gt;new($nb*2*2);
  my $i = $balancedhash-&gt;{$u};
  my $j = $balancedhash-&gt;{$v};
  $n-&gt;readNCode($net, $i, $j, $vcode);
  my $code = &amp;convertCodes($vcode);
  print STDERR "$u $i $v $j $code ". $vcode-&gt;to_Bin() . "\n";
  if ($code &gt; 0) {
    if (!defined $relations-&gt;{$code}) {
      $relations-&gt;{$code} = {};
    }
    $relations-&gt;{$code}-&gt;{$balanced-&gt;[$j]} = 1;
  }
  my $code = $n-&gt;readCode_1_1($i, $j);
  print STDERR "$u $i $v $j $code\n";
  return $relations;
}

sub getNCodesPair2 {
  my ($n, $net, $u, $v, $balanced, $balancedhash) = @_;
  my $num = $n-&gt;getNum();
  my $relations = {};
  my $nb = $n-&gt;{'numBits'};
  my $vcode = Bit::Vector-&gt;new($nb*2*2);
  my $i = $balancedhash-&gt;{$u};
  my $j = $balancedhash-&gt;{$v};
  $n-&gt;readNCode($net, $i, $j, $vcode);
  my $code = &amp;convertCodes($vcode);
  return $code;
}

sub getRelationship {
  my $counts = shift;
  if (scalar(grep { $counts-&gt;[$_]-&gt;[0] eq 0 } 1 .. 6) == 6) {
    return 0;
  }
  my $ratio = [];
  foreach my $c (0 .. 6) {
    if ($counts-&gt;[$c]-&gt;[1] &gt; 0) {
      $ratio-&gt;[$c] = $counts-&gt;[$c]-&gt;[0]/$counts-&gt;[$c]-&gt;[1];
    }
    else {
      $ratio-&gt;[$c] = 0;
    }
  }
  my $m1 = &amp;U::max([ map { $counts-&gt;[$_]-&gt;[0] } (2, 3, 5)]);
  my $m2 = &amp;U::max([ map { $counts-&gt;[$_]-&gt;[0] } (1, 4, 6)]);
  if ($m1 == $m2) {
    return 0;
  }
  my $m1 = &amp;U::max([ map { $ratio-&gt;[$_] } (2, 3, 5)]);
  my $m2 = &amp;U::max([ map { $ratio-&gt;[$_] } (1, 4, 6)]);
  if ($m1 == $m2) {
    return 0;
  }
  if ($m1 &gt; $m2) {
    if ($m1 == $ratio-&gt;[5] &amp;&amp; $m1 &gt; 0.5) {
      return 5;
    }
    if ($m1 == $ratio-&gt;[3] &amp;&amp; $m1 &gt; 0.5) {
      return 3;
    }
    if ($m1 == $ratio-&gt;[2] &amp;&amp; $m1 &gt; 0.5) {
      return 2;
    }
  }
  else {
    if ($m2 == $ratio-&gt;[6] &amp;&amp; $m2 &gt; 0.5) {
      return 6;
    }
    if ($m2 == $ratio-&gt;[4] &amp;&amp; $m2 &gt; 0.5) {
      return 4;
    }
    if ($m2 == $ratio-&gt;[1] &amp;&amp; $m2 &gt; 0.5) {
      return 1;
    }
  }
  return 0;
}

sub getRelationship2 {
  my ($counts, $i) = shift;
  if (scalar(grep { $counts-&gt;[$_]-&gt;[$i] eq 0 } 1 .. 6) == 6) {
    return 0;
  }
  my $ratio = [];
  foreach my $c (0 .. 6) {
    if ($counts-&gt;[$c]-&gt;[$i + 1] &gt; 0) {
      $ratio-&gt;[$c] = $counts-&gt;[$c]-&gt;[$i]/$counts-&gt;[$c]-&gt;[$i + 1];
    }
    else {
      $ratio-&gt;[$c] = 0;
    }
  }
  my $m1 = &amp;U::max([ map { $counts-&gt;[$_]-&gt;[$i] } (2, 3, 5)]);
  my $m2 = &amp;U::max([ map { $counts-&gt;[$_]-&gt;[$i] } (1, 4, 6)]);
  if ($m1 == $m2) {
    return 0;
  }
  my $m1 = &amp;U::max([ map { $ratio-&gt;[$_] } (2, 3, 5)]);
  my $m2 = &amp;U::max([ map { $ratio-&gt;[$_] } (1, 4, 6)]);
  if ($m1 == $m2) {
    return 0;
  }
  if ($m1 &gt; $m2) {
    if ($m1 == $ratio-&gt;[5] &amp;&amp; $m1 &gt; 0.5) {
      return 5;
    }
    if ($m1 == $ratio-&gt;[3] &amp;&amp; $m1 &gt; 0.5) {
      return 3;
    }
    if ($m1 == $ratio-&gt;[2] &amp;&amp; $m1 &gt; 0.5) {
      return 2;
    }
  }
  else {
    if ($m2 == $ratio-&gt;[6] &amp;&amp; $m2 &gt; 0.5) {
      return 6;
    }
    if ($m2 == $ratio-&gt;[4] &amp;&amp; $m2 &gt; 0.5) {
      return 4;
    }
    if ($m2 == $ratio-&gt;[1] &amp;&amp; $m2 &gt; 0.5) {
      return 1;
    }
  }
  return 0;
}

sub getRelationship3 {
  my $counts = shift;
  my $c1 = &amp;getRelationship2($counts, 0);
  my $c2 = &amp;getRelationship2($counts, 2);
  if ($c1 == $c2) {
    return $c1;
  }
  my $l = &amp;U::intersection([$c1, $c2], [2, 3, 5]);
  if (scalar(@$l) == 2) { return 5; }
  my $l = &amp;U::intersection([$c1, $c2], [1, 4, 6]);
  if (scalar(@$l) == 2) { return 6; }
  return 0;
}

sub getRel2 {
  my ($n, $net, $lu, $lv, $r1, $r2, $balanced, $balancedhash) = @_;
  my $nb = $n-&gt;{'numBits'};
  my $vcode = Bit::Vector-&gt;new($nb*2*2);
  foreach my $lid (@$lu) {
    my $li = $balancedhash-&gt;{$lid};
    foreach my $ljd (@$lv) {
      my $lj = $balancedhash-&gt;{$ljd};
      $n-&gt;readNCode($net, $li, $lj, $vcode);
      my $code = &amp;convertCodes($vcode);
      if ($code &gt; 0) {
        if (!defined $r1-&gt;{$code}) {
          $r1-&gt;{$code} = {};
        }
        if (!defined $r2-&gt;{$code}) {
          $r2-&gt;{$code} = {};
        }
        $r1-&gt;{$code}-&gt;{$balanced-&gt;[$li]} = 1;
        $r2-&gt;{$code}-&gt;{$balanced-&gt;[$lj]} = 1;
      }
    }
  }
}

sub getBooleanPath {
  my ($g, $nodes, $edges, $start, $rank, $ehash, $pathway, $nodelist) = @_;
  my $path = &amp;Graph::getDeepPath($nodes, $edges, $start, $rank);
  my $index = 0;
  my $ihash = {};
  my $ids = [];
  use UnionFind;
  my $uf = UnionFind-&gt;new();
  foreach my $id (map { $g-&gt;getNode($_) } @$path) {
    last if (defined $pathway-&gt;{$id});
    $pathway-&gt;{$id} = 1;
    my @k = keys %{$ehash-&gt;{$id}};
    foreach my $i (@k) {
      push @{$ihash-&gt;{$i}}, $index;
    }
    push @$ids, $id;
    $uf-&gt;makeset($id);
    $index++;
  }
  foreach my $i (keys %{$ihash}) {
    my $s = &amp;U::min($ihash-&gt;{$i});
    my $e = &amp;U::max($ihash-&gt;{$i});
    for (my $j = $s + 1; $j &lt;= $e; $j ++) {
      $uf-&gt;union($ids-&gt;[$j], $ids-&gt;[$j-1]);
    }
  }
  my $khash = {};
  foreach my $i (keys %{$ihash}) {
    my $s = &amp;U::min($ihash-&gt;{$i});
    my $e = &amp;U::max($ihash-&gt;{$i});
    for (my $j = $s; $j &lt;= $e; $j ++) {
      my $p = $uf-&gt;findset($ids-&gt;[$j]);
      $khash-&gt;{$p}-&gt;{$i} = 1;
    }
  }
  #print Dumper($ids);
  #print Dumper($ihash);
  #print Dumper($khash);
  my $clusters = {};
  my $bpath = [];
  foreach my $n (@$ids) {
    my $p = $uf-&gt;findset($n);
    if (!defined $clusters-&gt;{$p}) {
      push @$bpath, $p;
    }
    push @{$clusters-&gt;{$p}}, [$n, scalar(@{$nodelist-&gt;{$n}})];
  }
  my $res = [map { [$clusters-&gt;{$_}, [keys %{$khash-&gt;{$_}}]] } @$bpath];
  #print Dumper($res);
  return $res;
}

sub cleanupBooleanPath {
  my ($g, $path, $nodelist, $ehash) = @_;
  my ($nodes, $edges) = $g-&gt;getNodesEdges("5");
  my $kh = {};
  my $res = [];
  my $ekh = {};
  foreach my $k (@$path) {
    my $arr2 = [];
    foreach my $k2 (@{$k-&gt;[1]}) {
      next if (defined $ekh-&gt;{$k2});
      $ekh-&gt;{$k2} = 1;
      push @$arr2, $k2;
    }
    my $nl = [];
    foreach my $k1 (@{$k-&gt;[0]}) {
      my $ki = $g-&gt;getIndex($k1-&gt;[0]);
      next if (defined $kh-&gt;{$ki});
      push @$nl, $k1;
      $kh-&gt;{$ki} = 1;
    }
    next if (scalar(@$nl) == 0);
    my $extra = [];
    foreach my $k1 (@$nl) {
      my $i = $g-&gt;getIndex($k1-&gt;[0]);
      foreach my $ki (keys %{$edges-&gt;{$i}}) {
        next if (defined $kh-&gt;{$ki});
        my $n = $g-&gt;getNode($ki);
        next if (!defined $nodelist-&gt;{$n});
        if (defined $ehash-&gt;{$n}) {
          foreach my $k2 (keys %{$ehash-&gt;{$n}}) {
            next if (defined $ekh-&gt;{$k2});
            $ekh-&gt;{$k2} = 1;
            push @$arr2, $k2;
          }
        }
        push @$extra, [$n, scalar(@{$nodelist-&gt;{$n}})];
        $kh-&gt;{$ki} = 1;
      }
    }
    push @$res, [ [@$nl, @$extra], $arr2 ];
  }
  return $res;
}

sub getArijsNet {
  my $pre = "data/arijs-2018-uc";
  my $h = Hegemon-&gt;new(expr =&gt; "$pre\-expr.txt",
    idx =&gt; "$pre\-idx.txt", thr =&gt; "$pre\-thr.txt",
    survival =&gt; "$pre\-survival.txt");
  my $start = $h-&gt;getStart();
  my $end = $h-&gt;getEnd();
  my $all = [ $start .. $end ];
  use Network;
  my $file = "results/ibd-network-2.rl";
  my $n = Network-&gt;new(file =&gt; $file, mode =&gt; "r");
  $n-&gt;readNetworkFile();
  my $balanced = $n-&gt;getBalanced();
  my $balancedhash = {};
  for(my $i = 0; $i &lt; scalar(@{$balanced}); $i++) {
    $balancedhash-&gt;{$balanced-&gt;[$i]} = $i;
  }
  my $num = $n-&gt;getNum();
  my $rhash = {};
  my $net = $n-&gt;readAll();
  my $hash = &amp;U::getHash("data/peters-arjis-mapping.txt");
  my $res = { h =&gt; $h, n =&gt; $n, net =&gt; $net, hash =&gt; $hash,
    all =&gt; $all, balanced =&gt; $balanced, balancedhash =&gt; $balancedhash};
  return $res;
}

sub getArijsCode {
  my ($obj, $id1, $id2) = @_;
  my $all = $obj-&gt;{'all'};
  my $hash = $obj-&gt;{'hash'};
  my $h = $obj-&gt;{'h'};
  my $n = $obj-&gt;{'n'};
  my $net = $obj-&gt;{'net'};
  my $balanced = $obj-&gt;{'balanced'};
  my $balancedhash = $obj-&gt;{'balancedhash'};
  return [0, 0] if (!defined $hash-&gt;{$id1});
  return [0, 0] if (!defined $hash-&gt;{$id2});
  my $u = $hash-&gt;{$id1}-&gt;[1];
  my $v = $hash-&gt;{$id2}-&gt;[1];
  my $ex_u = $h-&gt;getExprData($u);
  my @exp_u = map { $ex_u-&gt;[$_] } @$all;
  my $ex_v = $h-&gt;getExprData($v);
  my @exp_v = map { $ex_v-&gt;[$_] } @$all;
  my $corr = &amp;U::correlation(\@exp_u, \@exp_v);
  my $code = &amp;getNCodesPair2($n, $net, $u, $v, $balanced, $balancedhash);
  return [$code, $corr]
}

sub isCodeCompatible {
  my ($code, $res) =  @_;
  if ($code &gt; 0) {
    return 1 if ($code eq $res-&gt;[0]);
    if ($res-&gt;[1] &gt; 0.5) {
      return 1 if ($code eq 2 || $code eq 3 || $code eq 5);
    }
    if ($res-&gt;[1] &lt; -0.25) {
      return 1 if ($code eq 1 || $code eq 4 || $code eq 6);
    }
  }
  return 0;
}

sub downloadData {
  my ($dbid, $pre) = @_;
  my $browser = LWP::UserAgent-&gt;new;
  my $url = "http://hegemon.ucsd.edu/Tools/explore.php";
  my $form = {'go' =&gt; 'dataDownload', 'id' =&gt; $dbid, 'genes' =&gt; '', 'groups' =&gt; '',
    'key' =&gt; 'polyps', 'param' =&gt; 'type:expr'};
  my $res = $browser-&gt;post($url, $form);
  my $ef = "$pre\-expr.txt";
  open(my $fh, "&gt;$ef");
  print $fh $res-&gt;content;
  close($fh);
  my $form = {'go' =&gt; 'dataDownload', 'id' =&gt; $dbid, 'genes' =&gt; '', 'groups' =&gt; '',
    'key' =&gt; 'polyps', 'param' =&gt; 'type:indexHeader'};
  my $res = $browser-&gt;post($url, $form);
  my $ef = "$pre\-ih.txt";
  open(my $fh, "&gt;$ef");
  print $fh $res-&gt;content;
  close($fh);
  $form-&gt;{'param'} = 'type:survival';
  my $res = $browser-&gt;post($url, $form);
  my $ef = "$pre\-survival.txt";
  open(my $fh, "&gt;$ef");
  print $fh $res-&gt;content;
  close($fh);
}

sub analyzeIBD {
  my @params = @_;
  if (defined $params[0] &amp;&amp; $params[0] eq "eq") {
    use Graph;
    my $g = Graph-&gt;new();
    my $file = "results/ibd-network-res.txt";
    my $count = 0;
    open(my $fh, "&lt;$file");
    while (&lt;$fh&gt;) {
      next if (!/^Found : 5/);
      chomp;
      my @list = split("\t");
      my $id1 = $list[3];
      my $id2 = $list[4];
      $g-&gt;addEdge($id1, $id2, "eq");
      $count++;
      #if ($count &gt; 10000) {
      #  last;
      #}
    }
    close($fh);
    my ($nodes, $edges) = $g-&gt;getNodesEdges("eq");
    my $hash = {};
    my $shash = {};
    use UnionFind;
    my $uf = UnionFind-&gt;new();
    my $count = 0;
    my $num = scalar(keys %{$edges});
    foreach my $u (keys %{$edges}) {
      print STDERR "$count/$num\n";
      my $data = [];
      my $scores = [];
      foreach my $v (keys %{$edges-&gt;{$u}}) {
        my $r = &amp;rho($u, $v, $edges, $hash, $shash);
        push @$data, $r;
        push @$scores, [$r, $v];
        $uf-&gt;makeset($u);
        $uf-&gt;makeset($v);
      }
      my $filter = [sort { $b-&gt;[0] &lt;=&gt; $a-&gt;[0] } @$scores];
      foreach my $s (@$filter) {
        if ($uf-&gt;findset($u) != $uf-&gt;findset($s-&gt;[1])) {
          $uf-&gt;union($u, $s-&gt;[1]);
          print join("\t", $g-&gt;getNode($u), $s-&gt;[0], 
            $g-&gt;getNode($s-&gt;[1])), "\n";
          last;
        }
      }
      $count++;
    }
  }
  if (defined $params[0] &amp;&amp; $params[0] eq "eq-corr") {
    my $pre = "data/arijs-2018-uc";
    my $h = Hegemon-&gt;new(expr =&gt; "$pre\-expr.txt",
      idx =&gt; "$pre\-idx.txt", thr =&gt; "$pre\-thr.txt",
      survival =&gt; "$pre\-survival.txt");
    my $start = $h-&gt;getStart();
    my $end = $h-&gt;getEnd();
    my $all = [ $start .. $end ];
    my $hash = &amp;U::getHash("data/peters-arjis-mapping.txt");
    my $file = "data/ibd-network-g-eq-4.txt";
    my $index = 0;
    open(my $fh, "&lt;$file");
    while (&lt;$fh&gt;) {
      chomp;
      my @list = split("\t");
      my $id1 = $list[0];
      my $id2 = $list[2];
      next if (!defined $hash-&gt;{$id1});
      next if (!defined $hash-&gt;{$id2});
      my $u = $hash-&gt;{$id1}-&gt;[1];
      my $v = $hash-&gt;{$id2}-&gt;[1];
      my $ex_u = $h-&gt;getExprData($u);
      my @exp_u = map { $ex_u-&gt;[$_] } @$all;
      my $ex_v = $h-&gt;getExprData($v);
      my @exp_v = map { $ex_v-&gt;[$_] } @$all;
      my $corr = &amp;U::correlation(\@exp_u, \@exp_v);
      print join("\t", $id1, $id2, $list[1], $corr), "\n";
      $index++;
    }
    close($fh);
  }
  if (defined $params[0] &amp;&amp; $params[0] eq "eq-corr-thr") {
    my $file = "results/ibd-network-g-eq-4-corr.txt";
    my $index = 0;
    open(my $fh, "&lt;$file");
    my $data = [];
    while (&lt;$fh&gt;) {
      chomp;
      my @list = split("\t");
      if ($list[3] &gt; 0) {
        push @$data, $list[3];
      }
      $index++;
    }
    close($fh);
    $data = [sort { $a &lt;=&gt; $b } @$data];
    my $res = &amp;U::fitStep($data, 0, scalar(@$data) - 1);
    my $thr = [$res-&gt;[6], $res-&gt;[3], $res-&gt;[6]-0.5, $res-&gt;[6]+0.5];
    print STDERR join("\t", "DR", @$thr), "\n";
  }
  if (defined $params[0] &amp;&amp; $params[0] eq "cls") {
    use UnionFind;
    my $uf = UnionFind-&gt;new();
    my $nodes = {};
    my $hash = &amp;U::getHash("data/peters-arjis-mapping.txt");
    use Network;
    my $file = "results/ibd-network-2.rl";
    my $n = Network-&gt;new(file =&gt; $file, mode =&gt; "r");
    $n-&gt;readNetworkFile();
    my $balanced = $n-&gt;getBalanced();
    my $balancedhash = {};
    for(my $i = 0; $i &lt; scalar(@{$balanced}); $i++) {
      $balancedhash-&gt;{$balanced-&gt;[$i]} = $i;
      my $id = $balanced-&gt;[$i];
      if (defined $hash-&gt;{$id}) {
        $uf-&gt;makeset($id);
        $nodes-&gt;{$id} = {} if (!defined $nodes-&gt;{$id});
        $nodes-&gt;{$id}-&gt;{$id} = 1;
      }
    }
    my $num = $n-&gt;getNum();
    my $rhash = {};
    my $net = $n-&gt;readAll();
    my $file = "data/ibd-network-g-eq-4.txt";
    open(my $fh, "&lt;$file");
    while (&lt;$fh&gt;) {
      chomp;
      my @list = split("\t");
      my $id1 = $list[0];
      my $id2 = $list[2];
      if (defined $hash-&gt;{$id1}) {
        $uf-&gt;makeset($id1);
        $nodes-&gt;{$id1} = {} if (!defined $nodes-&gt;{$id1});
        $nodes-&gt;{$id1}-&gt;{$id1} = 1;
      }
      if (defined $hash-&gt;{$id2}) {
        $uf-&gt;makeset($id2);
        $nodes-&gt;{$id2} = {} if (!defined $nodes-&gt;{$id2});
        $nodes-&gt;{$id2}-&gt;{$id2} = 1;
      }
    }
    close($fh);
    my $file = "results/ibd-network-g-eq-4-corr.txt";
    my $index = 0;
    open(my $fh, "&lt;$file");
    while (&lt;$fh&gt;) {
      chomp;
      my @list = split("\t");
      my $id1 = $list[0];
      my $id2 = $list[1];
      next if (!defined $hash-&gt;{$id1});
      next if (!defined $hash-&gt;{$id2});
      $uf-&gt;makeset($id1);
      $uf-&gt;makeset($id2);
      $nodes-&gt;{$id1} = {} if (!defined $nodes-&gt;{$id1});
      $nodes-&gt;{$id2} = {} if (!defined $nodes-&gt;{$id2});
      $nodes-&gt;{$id1}-&gt;{$id1} = 1;
      $nodes-&gt;{$id2}-&gt;{$id2} = 1;
      my $u = $hash-&gt;{$id1}-&gt;[1];
      my $v = $hash-&gt;{$id2}-&gt;[1];
      if ($list[3] &gt; 0.5) {
        push @$data, $list[3];
        $uf-&gt;union($id1, $id2);
        $nodes-&gt;{$id1}-&gt;{$id2} = $list[1];
        $nodes-&gt;{$id2}-&gt;{$id1} = $list[1];
      }
      elsif (defined $balancedhash-&gt;{$u} &amp;&amp; defined $balancedhash-&gt;{$v}) {
        my $code = &amp;getNCodesPair2($n, $net, $u, $v, $balanced, $balancedhash);
        next if ($code ne 5 &amp;&amp; $code ne 3 &amp;&amp; $code ne 2);
        push @$data, $list[3];
        $uf-&gt;union($id1, $id2);
        $nodes-&gt;{$id1}-&gt;{$id2} = $list[3];
        $nodes-&gt;{$id2}-&gt;{$id1} = $list[3];
      }
      $index++;
    }
    close($fh);
    my $rank = {};
    foreach my $n (keys %{$nodes}) {
      $rank-&gt;{$n} = scalar(keys %{$nodes-&gt;{$n}});
    }
    my $clusters = {};
    foreach my $n (keys %{$nodes}) {
      push @{$clusters-&gt;{$uf-&gt;findset($n)}}, $n;
    }
    if (defined $params[1]) {
      my $n = $uf-&gt;findset($params[1]);
      my @l = sort { $rank-&gt;{$b} &lt;=&gt; $rank-&gt;{$a} } @{$clusters-&gt;{$n}};
      print scalar(@l), "\n";
      foreach my $id (@l) {
        print join("\t", $id, $rank-&gt;{$id}), "\n";
      }
    }
    else {
      foreach my $n (keys %{$clusters}) {
        my @l = sort { $rank-&gt;{$b} &lt;=&gt; $rank-&gt;{$a} } @{$clusters-&gt;{$n}};
        print join("\t", $l[0], scalar(@l), @l), "\n";
      }
    }
  }
  if (defined $params[0] &amp;&amp; $params[0] eq "g-cls") {
    my $file = "results/ibd-network-g-eq-cls-4.txt";
    my $nodes = {};
    my $ids = [];
    open(my $fh, "&lt;$file");
    while (&lt;$fh&gt;) {
      chomp;
      my @list = split("\t");
      my $id = $list[0];
      my $l = [@list[2 .. $#list]];
      $nodes-&gt;{$id} = $l;
      push @$ids, $id;
    }
    close($fh);
    use Network;
    my $file = "results/ibd-network-1.rl";
    my $n = Network-&gt;new(file =&gt; $file, mode =&gt; "r");
    $n-&gt;readNetworkFile();
    my $balanced = $n-&gt;getBalanced();
    my $balancedhash = {};
    for(my $i = 0; $i &lt; scalar(@{$balanced}); $i++) {
      $balancedhash-&gt;{$balanced-&gt;[$i]} = $i;
    }
    my $num = $n-&gt;getNum();
    my $rhash = {};
    my $net = $n-&gt;readAll();
    my $arijs = &amp;getArijsNet();
    foreach my $i (0 .. $#{$ids}) {
      my $u = $ids-&gt;[$i];
      if (defined $params[1]) {
        $u = $params[1];
      }
      print STDERR "$i $u\n";
      my $n1 = scalar(@{$nodes-&gt;{$u}});
      my $mid = int($n1 / 2);
      my $l = [$u];
      if ($n1 &gt; 2) {
        $l = [ $u, map { $nodes-&gt;{$u}-&gt;[$_] } (1, $mid)];
      }
      if ($n1 &gt; 10) {
        push @$l, map { $nodes-&gt;{$u}-&gt;[$_] } (int($mid/4), int($mid/2), $mid - 1);
      }
      my $ru = &amp;getNCodesArijs($arijs, $n, $net, $l, $balanced, $balancedhash);
      if (defined $params[2]) {
        my $v = $params[2];
        print STDERR $v, "\n";
        foreach my $u (@$l) {
          print STDERR $u, "\n";
          &amp;getNCodesPair($n, $net, $u, $v, $balanced, $balancedhash);
        }
      }
      foreach my $c (1 .. 6) {
        foreach my $v (keys %{$ru-&gt;{$c}}) {
          next if (!defined $nodes-&gt;{$v});
          next if ($n1 &gt; 10 &amp;&amp; $ru-&gt;{$c}-&gt;{$v} &lt; 3);
          next if (defined $params[2] &amp;&amp; $params[2] ne $v);
          print join("\t", $u, $c, $v, $ru-&gt;{$c}-&gt;{$v}), "\n";
        }
      }
      if (defined $params[1]) {
        last;
      }
    }
  }
  if (defined $params[0] &amp;&amp; $params[0] eq "g-filter") {
    my $arijs = &amp;getArijsNet();
    my $file = "results/ibd-network-g-eq-g-4.txt";
    my $index = 0;
    open(my $fh, "&lt;$file");
    while (&lt;$fh&gt;) {
      chomp;
      my @list = split("\t");
      my $id1 = $list[0];
      my $id2 = $list[2];
      my $res = &amp;getArijsCode($arijs, $id1, $id2);
      if (&amp;isCodeCompatible($list[1], $res)) {
        print join("\t", @list, $res-&gt;[1]), "\n";
      }
      $index++;
      if ( ($index % 10000) == 0 ) {
        print STDERR "$index\n";
      }
    }
    close($fh);
  }
  if (defined $params[0] &amp;&amp; $params[0] eq "genes") {
    use Graph;
    my $g = Graph-&gt;new();
    my $file = "data/ibd-network-g-eq-g-4.txt";
    my $nodelist = {};
    my $ids = [];
    open(my $fh, "&lt;$file");
    while (&lt;$fh&gt;) {
      chomp;
      my @list = split("\t");
      $g-&gt;addEdge($list[0], $list[2], $list[1]);
    }
    close($fh);
    my $file = "data/ibd-network-g-eq-cls-4.txt";
    my $nhash = {};
    my $rank = {};
    my $ids = [];
    open(my $fh, "&lt;$file");
    while (&lt;$fh&gt;) {
      chomp;
      my @list = split("\t");
      my $id = $list[0];
      my $l = [@list[2 .. $#list]];
      $nodelist-&gt;{$id} = $l;
      $nhash-&gt;{$id} = $id;
      foreach my $i (@$l) {
        $nhash-&gt;{$i} = $id;
      }
      $rank-&gt;{$g-&gt;getIndex($id)} = scalar(@$l);
      push @$ids, $g-&gt;getIndex($id);
    }
    close($fh);
    my $ehash = {};
    my ($nodes, $edges) = $g-&gt;getNodesEdges("5");
    if (defined $params[1]) {
      foreach my $id (@params) {
        next if (!defined $nhash-&gt;{$id});
        my $i = $g-&gt;getIndex($nhash-&gt;{$id});
        $ehash-&gt;{$g-&gt;getNode($i)}-&gt;{$id} = 1;
        if ($rank-&gt;{$i} &lt; 2000) {
          $rank-&gt;{$i} += 5000;
        }
        foreach my $k (keys %{$edges-&gt;{$i}}) {
          if ($rank-&gt;{$k} &gt; 5 &amp;&amp; $rank-&gt;{$k} &lt; 2000) {
            $rank-&gt;{$k} += 2000;
          }
          $ehash-&gt;{$g-&gt;getNode($k)}-&gt;{$id} = 1;
        }
      }
    }
    #print Dumper($ehash);
    $ids = [sort { $rank-&gt;{$b} &lt;=&gt;  $rank-&gt;{$a} } @$ids];
    my $res = [];
    push @$res, [0, scalar(@$ids), 
      [map { [$_, $rank-&gt;{$g-&gt;getIndex($_)}] } 
      map { $g-&gt;getNode($ids-&gt;[$_]) } 0 .. 5], []];
    my $pathway = {};
    my ($nodes, $edges) = $g-&gt;getNodesEdges("2");
    my $start = $ids-&gt;[0];
    my $path = &amp;getBooleanPath($g, $nodes, $edges, $start, $rank, $ehash, $pathway, $nodelist);
    my $path = &amp;cleanupBooleanPath($g, $path, $nodelist, $ehash);
    my $index = 1;
    foreach my $k (@$path) {
      push @$res, [-$index, scalar(@{$k-&gt;[0]}), @$k];
      $index++;
    }
    my ($nodes, $edges6) = $g-&gt;getNodesEdges("6");
    my ($nodes, $edges4) = $g-&gt;getNodesEdges("4");
    my $edges = &amp;Graph::doUnion($edges6, $edges4);
    my $ids = [sort { $rank-&gt;{$b} &lt;=&gt; $rank-&gt;{$a} } keys %{$edges-&gt;{$start}}];
    push @$res, [0, scalar(@$ids), 
      [map { [$_, $rank-&gt;{$g-&gt;getIndex($_)}] } 
      map { $g-&gt;getNode($ids-&gt;[$_]) } 0 .. 5], []];
    my ($nodes, $edges) = $g-&gt;getNodesEdges("2");
    my $start = $ids-&gt;[0];
    my $path = &amp;getBooleanPath($g, $nodes, $edges, $start, $rank, $ehash, $pathway, $nodelist);
    my $path = &amp;cleanupBooleanPath($g, $path, $nodelist, $ehash);
    my $index = 1;
    foreach my $k (@$path) {
      push @$res, [$index, scalar(@{$k-&gt;[0]}), @$k];
      $index++;
    }

    my $json = encode_json $res;
    print $json, "\n";
    #foreach (@$res) {
    #  print join("\t", @$_), "\n";
    #}
  }
  if (defined $params[0] &amp;&amp; $params[0] eq "query") {
    my $file = "results/ibd-network-g-eq-cls-4.txt";
    my $nodelist = {};
    my $nhash = {};
    my $rank = {};
    my $ids = [];
    open(my $fh, "&lt;$file");
    while (&lt;$fh&gt;) {
      chomp;
      my @list = split("\t");
      my $id = $list[0];
      my $l = [@list[2 .. $#list]];
      $nodelist-&gt;{$id} = $l;
      $nhash-&gt;{$id} = $id;
      foreach my $i (@$l) {
        $nhash-&gt;{$i} = $id;
      }
      $rank-&gt;{$id} = scalar(@$l);
      push @$ids, $id;
    }
    close($fh);
    use Network;
    my $file = "results/ibd-network-1.rl";
    my $n = Network-&gt;new(file =&gt; $file, mode =&gt; "r");
    $n-&gt;readNetworkFile();
    my $balanced = $n-&gt;getBalanced();
    my $balancedhash = {};
    for(my $i = 0; $i &lt; scalar(@{$balanced}); $i++) {
      $balancedhash-&gt;{$balanced-&gt;[$i]} = $i;
    }
    my $num = $n-&gt;getNum();
    my $rhash = {};
    my $net = $n-&gt;readAll();
    my $u = $nhash-&gt;{$params[1]};
    my $n1 = scalar(@{$nodelist-&gt;{$u}});
    print "$u ($n1)\n";
    my $l = $nodelist-&gt;{$u};
    my $v = $nhash-&gt;{$params[2]};
    my $n2 = scalar(@{$nodelist-&gt;{$v}});
    print "$v ($n2)\n";
    my $ru = &amp;getNCodes($n, $net, $l, $balanced, $balancedhash);
    my $vhash = {};
    foreach my $v1 (@{$nodelist-&gt;{$v}}) {
      $vhash-&gt;{$v1} = 1;
    }
    my @rels = (1 .. 6);
    if (defined $params[3]) {
      @rels = grep { $_ =~ /[0-9]/ } split("", $params[3]);
    }
    foreach my $c (@rels) {
      foreach my $v (keys %{$ru-&gt;{$c}}) {
        next if (!defined $vhash-&gt;{$v});
        print join("\t", $u, $c, $v, $ru-&gt;{$c}-&gt;{$v}), "\n";
      }
    }
  }
  if (defined $params[0] &amp;&amp; $params[0] eq "path") {
    use Graph;
    my $g = Graph-&gt;new();
    my $file = "results/ibd-network-g-eq-g-4.txt";
    my $nodelist = {};
    my $ids = [];
    open(my $fh, "&lt;$file");
    while (&lt;$fh&gt;) {
      chomp;
      my @list = split("\t");
      $g-&gt;addEdge($list[0], $list[2], $list[1]);
    }
    close($fh);
    my $file = "results/ibd-network-g-eq-cls-4.txt";
    my $nhash = {};
    my $rank = {};
    my $ids = [];
    open(my $fh, "&lt;$file");
    while (&lt;$fh&gt;) {
      chomp;
      my @list = split("\t");
      my $id = $list[0];
      my $l = [@list[2 .. $#list]];
      $nodelist-&gt;{$id} = $l;
      $nhash-&gt;{$id} = $id;
      foreach my $i (@$l) {
        $nhash-&gt;{$i} = $id;
      }
      $rank-&gt;{$g-&gt;getIndex($id)} = scalar(@$l);
      push @$ids, $g-&gt;getIndex($id);
    }
    close($fh);
    my $u = $nhash-&gt;{$params[1]};
    my $n1 = scalar(@{$nodelist-&gt;{$u}});
    print "$u ($n1)\n";
    my $l = $nodelist-&gt;{$u};
    my $v = $nhash-&gt;{$params[2]};
    my $n2 = scalar(@{$nodelist-&gt;{$v}});
    print "$v ($n2)\n";
    my $res = $g-&gt;Path("5", $u, $v, 0, 100);
    print join("\t", map { $g-&gt;getNode($res-&gt;[$_])." ".
     $g-&gt;getType($res-&gt;[$_], $res-&gt;[$_ + 1]) } 0 .. $#{$res}), "\n";
  }
  if (defined $params[0] &amp;&amp; $params[0] eq "rank") {
    my $file = "results/ibd-network-g-eq-cls-4.txt";
    my $nodelist = {};
    my $nhash = {};
    my $rank = {};
    my $ids = [];
    open(my $fh, "&lt;$file");
    while (&lt;$fh&gt;) {
      chomp;
      my @list = split("\t");
      my $id = $list[0];
      my $l = [@list[2 .. $#list]];
      $nodelist-&gt;{$id} = $l;
      $nhash-&gt;{$id} = $id;
      foreach my $i (@$l) {
        $nhash-&gt;{$i} = $id;
      }
      $rank-&gt;{$id} = scalar(@$l);
      push @$ids, $id;
    }
    close($fh);
    my $u = $nhash-&gt;{$params[1]};
    my $n1 = scalar(@{$nodelist-&gt;{$u}});
    print "$u ($n1)\n";
    my $l = $nodelist-&gt;{$u};
    my $v = $nhash-&gt;{$params[2]};
    my $n2 = scalar(@{$nodelist-&gt;{$v}});
    print "$v ($n2)\n";
    my @rels = (1 .. 6);
    if (defined $params[3]) {
      @rels = grep { $_ =~ /[0-9]/ } split("", $params[3]);
    }
    my $vhash = {};
    if (defined $params[4]) {
      $vhash = &amp;U::getHash($params[4], 1);
    }
    else {
      foreach my $v1 (@{$nodelist-&gt;{$v}}) {
        $vhash-&gt;{$v1} = 1;
      }
    }
    use Network;
    my $file = "results/ibd-network-1.rl";
    my $n = Network-&gt;new(file =&gt; $file, mode =&gt; "r");
    $n-&gt;readNetworkFile();
    my $balanced = $n-&gt;getBalanced();
    my $balancedhash = {};
    for(my $i = 0; $i &lt; scalar(@{$balanced}); $i++) {
      $balancedhash-&gt;{$balanced-&gt;[$i]} = $i;
    }
    my $num = $n-&gt;getNum();
    my $rhash = {};
    my $net = $n-&gt;readAll();
    my $ru = &amp;getNCodes($n, $net, $l, $balanced, $balancedhash);
    foreach my $c (@rels) {
      foreach my $v (keys %{$ru-&gt;{$c}}) {
        next if (!defined $vhash-&gt;{$v});
        print join("\t", $u, $c, $v, $ru-&gt;{$c}-&gt;{$v}), "\n";
      }
    }
  }
}

</pre></body></html>