计算图中子图的出现次数

Question

我有一个 3 维数据集，描述了可以用图表表示的基因相互作用。数据集样本为：

a + b  
b + c  
c - f  
b - d  
a + c  
f + g  
g + h  
f + h  

“+”表示左侧基因正向调节右侧基因。在这些数据中，我想计算一个基因（例如，x）正向调节另一个基因（例如，y），y 反过来正向调节另一个基因（例如，z）的子图。此外，z也受到x的正向调节。上图中有两种这样的情况。我想最好使用 awk 来执行此搜索，但任何脚本语言都可以。我很抱歉提出了一个过于具体的问题，并提前感谢您的帮助。

Answer 1

注意：请参阅下面有关 Graphviz 的信息。

这应该为您提供一个起点：

编辑：此版本处理由多个字符描述的基因。

awk '
    BEGIN { regdelim = "|" }
    {
        delim=""
        if ($2 == "+") {
            if (plus[$1]) delim=regdelim
            plus[$1]=plus[$1] delim $3
        }
        else
            if ($2 == "-") {
            if (minus[$1]) delim=regdelim
                minus[$1]=minus[$1] delim $3
            }
    }
    END {
        for (root in plus) {
            split(plus[root],regs,regdelim)
            for (reg in regs) {
                if (plus[regs[reg]] && plus[root] ~ plus[regs[reg]]) {
                    print "Match: ", root, "+", regs[reg], "+", plus[regs[reg]]
                }
            }
        }
    }
' inputfile

在 BEGIN 子句中，将 regdelim 设置为数据中未出现的字符。

我省略了负数据的处理代码。

输出：

Match:  a + b + c
Match:  f + g + h

编辑 2：

下面的版本允许您搜索任意组合。它概括了原始版本中使用的技术，因此不需要重复代码。它还修复了其他一些错误限制。

#!/bin/bash
# written by Dennis Williamson - 2010-11-12
# for http://stackoverflow.com/questions/4161001/counting-the-occurrence-of-a-sub-graph-in-a-graph
# A (AB) B, A (AC) C, B (BC) C - where "(XY)" represents a + or a - 
# provided by the positional parameters $1, $2 and $3
# $4 carries the data file name and is referenced at the end of the script
awk -v AB=$1 -v AC=$2 -v BC=$3 '
    BEGIN { regdelim = "|" }
    {
        if ($2 == AB) {
            if (regAB[$1]) delim=regdelim; else delim=""
            regAB[$1]=regAB[$1] delim $3
        }
        if ($2 == AC) {
            if (regAC[$1]) delim=regdelim; else delim=""
            regAC[$1]=regAC[$1] delim $3
        }
        if ($2 == BC) {
            if (regBC[$1]) delim=regdelim; else delim=""
            regBC[$1]=regBC[$1] delim $3
        }
    }
    END {
        for (root in regAB) {
            split(regAB[root],ABarray,regdelim)
            for (ABindex in ABarray) {
                split(regAC[root],ACarray,regdelim)
                for (ACindex in ACarray) {
                    split(regBC[ABarray[ABindex]],BCarray,regdelim)
                    for (BCindex in BCarray) {
                        if (ACarray[ACindex] == BCarray[BCindex]) {
                            print "    Match:", root, AB, ABarray[ABindex] ",", root, AC, ACarray[ACindex] ",", ABarray[ABindex], BC, BCarray[BCindex]
                        }
                    }
                }
            }
        }
    }
' "$4"

可以这样调用它来进行详尽的搜索：

for ab in + -; do for ac in + -; do for bc in + -; do echo "Searching: $ab$ac$bc"; ./searchgraph $ab $ac $bc inputfile; done; done; done

对于此数据：

a - e
a + b
b + c
c - f
m - n
b - d
a + c
b - e
l - n
f + g
b + i
g + h
l + m
f + h
a + i
a - j
k - j
a - k

调用新版本脚本的 shell 循环的输出将如下所示：

Searching: +++
    Match: a + b, a + c, b + c
    Match: a + b, a + i, b + i
    Match: f + g, f + h, g + h
Searching: ++-
Searching: +-+
Searching: +--
    Match: l + m, l - n, m - n
    Match: a + b, a - e, b - e
Searching: -++
Searching: -+-
Searching: --+
Searching: ---
    Match: a - k, a - j, k - j

Edit 3:

Graphviz

另一种方法是使用 Graphviz。 DOT 语言可以描述图形和gvpr，这是一种“类似 AWK”¹ 编程语言，可以分析和操作 DOT 文件。

给定问题中所示格式的输入数据，您可以使用以下 AWK 程序将其转换为 DOT：

#!/usr/bin/awk -f
BEGIN {
    print "digraph G {"
    print "    size=\"5,5\""
    print "    ratio=.85"
    print "    node [fontsize=24 color=blue penwidth=3]"
    print "    edge [fontsize=18 labeldistance=5 labelangle=-8 minlen=2 penwidth=3]"
    print "    {rank=same; f l}"
    m  = "-"    # ASCII minus/hyphen as in the source data
    um = "−"    # u2212 minus: − which looks better on the output graphic
    p  = "+"
}

{
    if ($2 == m) { $2 = um; c = lbf = "red"; arr=" arrowhead = empty" }
    if ($2 == p) { c = lbf = "green3"; arr="" }
    print "    " $1, "->", $3, "[taillabel = \"" $2 "\" color = \"" c "\" labelfontcolor = \"" lbf "\"" arr "]"
}
END {
    print "}"
}

要运行的命令如下所示：

$ ./dat2dot data.dat > data.dot

然后您可以使用以下命令创建上面的图形：

$ dot -Tpng -o data.png data.dot

我在这个答案中使用了上面给出的扩展数据。

要对指定的子图类型进行详尽的搜索，您可以使用以下 gvpr 程序：

BEGIN {
    edge_t AB, BC, AC;
}

E {
    AB = $;
    BC = fstedge(AB.head);
    while (BC && BC.head.name != AB.head.name) {
        AC = isEdge(AB.tail,BC.head,"");
        if (AC) {
            printf("%s %s %s, ", AB.tail.name, AB.taillabel, AB.head.name);
            printf("%s %s %s, ", AC.tail.name, AC.taillabel, AC.head.name);
            printf("%s %s %s\n", BC.tail.name, BC.taillabel, BC.head.name);
        }
        BC = nxtedge(BC, AB.head);
    }
}

要运行它，您可以使用：

$ gvpr -f groups.g data.dot | sort -k 2,2 -k 5,5 -k 8,8

输出将与上面的 AWK/shell 组合的输出类似（在“编辑 2”下）：

a + b, a + c, b + c
a + b, a + i, b + i
f + g, f + h, g + h
a + b, a − e, b − e
l + m, l − n, m − n
a − k, a − j, k − j

---

¹ _{宽松地说。}

Answer 2

下面是使用 Perl 的非常规方法。

#! /usr/bin/perl

use warnings;
use strict;

my $graph = q{
  a + c
  b + c
  c - f
  b - d
  a + b
  f + g
  g + h
  f + h
};

my $nodes = join ",", sort keys %{ { map +($_ => 1), $graph =~ /(\w+)/g } };
my $search = "$nodes:$nodes:$nodes:$graph";

my $subgraph = qr/
  \A  .*?  (?<x>\w+)  .*:
      .*?  (?<y>\w+)  .*:
      .*?  (?<z>\w+)  .*:
  (?= .*^\s*  \k<x>  \s*  \+  \s*  \k<y>  \s*$)
  (?= .*^\s*  \k<y>  \s*  \+  \s*  \k<z>  \s*$)
  (?= .*^\s*  \k<x>  \s*  \+  \s*  \k<z>  \s*$)
  (?{ print "x=$+{x}, y=$+{y}, z=$+{z}\n" })
  (?!)
/smx;

$search =~ /$subgraph/;

正则表达式引擎是一个强大的工具。对于您的问题，我们描述传递子图的结构，然后允许生成的机器查找所有子图。

输出：

x=a, y=b, z=c
x=f, y=g, z=h

---

使用相同技术的更通用的工具是可能的。例如，假设您希望能够指定基因模式，例如 a+b+c;a+c 或 g1+g2-g3;g1+g3。我希望这些模式的含义是显而易见的。

在前面的内容中，我指定了最低版本 5.10.0，因为代码使用 // 和词法 $_。该代码构造了将评估代码的动态正则表达式，这是 use re 'eval' pragma 启用的。

#! /usr/bin/perl

use warnings;
use strict;

use 5.10.0;
use re 'eval';

标识符是一个或多个“单词字符”的序列，即字母、数字或下划线。

my $ID = qr/\w+/;

给定一个接受变量名称的正则表达式，unique_vars 会在某个规范中搜索所有变量名称并返回它们而不重复。

sub unique_vars {
  my($_,$pattern) = @_;
  keys %{ { map +($_ => undef), /($pattern)/g } };
}

将基因模式编译成正则表达式有点麻烦。它动态生成一个搜索目标和正则表达式，其形式与上面的静态形式相同。

多次出现逗号分隔变量的第一部分让正则表达式引擎尝试每个基因的每个可能值。然后，前瞻 (?=...) 扫描图形，查找具有所需属性的边。如果所有前瞻都成功，我们就会记录命中。

最后的奇怪的正则表达式 (?!) 是一个无条件失败，迫使匹配器回溯并尝试与不同的基因进行匹配。因为它是无条件的，所以引擎将评估所有可能性。

同时从多个线程调用同一个闭包可能会产生奇怪的结果。

sub compile_gene_pattern {
  my($dataset,$pattern) = @_;
  my @vars   = sort +unique_vars $pattern, qr/[a-z]\d*/;  # / for SO hilite
  my $nodes  = join ",", sort +unique_vars $dataset, $ID;
  my $search = join("", map "$_:", ($nodes) x @vars) . "\n"
             . $dataset;

  my $spec = '\A' . "\n" . join("", map ".*?  (?<$_>$ID)  .*:\n", @vars);
  for (split /;/, $pattern) {
    while (s/^($ID)([-+])($ID)/$3/) {
      $spec .= '(?= .*^\s*  ' .
               ' \b\k<' .           $1  . '>\b ' .
               ' \s*'   . quotemeta($2) . '\s* ' .
               ' \b\k<' .           $3  . '>\b ' .
               ' \s*$)' . "\n";
    }
  }
  my %hits;
  $spec .= '(?{ ++$hits{join "-", @+{@vars}} })' . "\n" .
           '(?!) # backtrack' . "\n";

  my $nfa = eval { qr/$spec/smx } || die "$0: INTERNAL: bad regex:\n$@";
  sub {
    %hits = ();  # thread-safety? :-(
    (my $_ = $search) =~ /$nfa/;
    map [split /-/], sort keys %hits;
  }
}

阅读数据集并让用户了解任何问题。

sub read_dataset {
  my($path) = @_;

  open my $fh, "<", $path or die "$0: open $path: $!";

  local $/ = "\n";
  local $_;
  my $graph;

  my @errors;
  while (<$fh>) {
    next if /^\s*#/ || /^\s*$/;

    if (/^ \s* $ID \s* [-+] \s* $ID \s* $/x) {
      $graph .= $_;
    }
    else {
      push @errors, "$.: syntax error";
    }
  }

  return $graph unless @errors;

  die map "$0: $path:$_\n", @errors;
}

现在我们将一切付诸实践：

my $graphs = shift // "graphs.txt";
my $dataset = read_dataset $graphs;

my $ppp = compile_gene_pattern $dataset, "a+b+c;a+c";
print "@$_\n" for $ppp->();

my $pmp = compile_gene_pattern $dataset, "g1+g2-g3;g1+g3";
print "@$_\n" for $pmp->();

给定 graphs.txt 的内容

a + b  
b + c  
c - f  
b - d  
a + c  
f + g  
g + h  
f + h

foo + bar
bar - baz
foo + baz

，然后运行程序，我们得到以下输出：

a b c
f g h
foo bar baz

Answer 3

我假设“计算子图”的意思是计算子图中的节点。如果这就是您需要的，您可以使用任何脚本语言，并且必须存储图形，首先，通过创建存储图形的结构或类，节点结构/类应该如下所示（这不符合任何语言的语法，这只是您的应用程序的计划）：

Node {color = 0; title = ""; minusNodeSet = null; plusNodeSet = null}

其中 color = 0（颜色的默认值意味着您之前没有访问过该节点），标题将为 'a', 'b', 'c '， 等等。 minusNodeSet 是存储这些节点的节点集，其中负顶点指向我们的节点，plusNodeSet 是存储这些节点的节点集，其中正顶点指向我们的节点。

现在，我们有了一个架构，应该在深度优先算法中使用它：

int depth_first(Node actualNode)
{
    if (actualNode.color == 1)
    return;

    number = 1;
    actualNode.color = 1;
    foreach actualNode.nodeSet as node do
        if (node.color == 0)
            number = number + depth_first(node);
    return number;
}

如果我误解了你的问题，请告诉我，以便能够将我的答案编辑为更有用的答案。

Answer 4

我的其他中的正则表达式的结构答案类似于列表单子处理。受此启发，对传递子图的搜索如下作为 literate Haskell。将此答案复制并粘贴到扩展名为 .lhs 的文件中以获取工作程序。请务必用空行将由前导 > 标记的代码部分括起来。

感谢您提出有趣的问题！

前面的一点：

> {-# LANGUAGE ViewPatterns #-}

> module Main where

> import Control.Monad (guard)
> import Data.List (nub)
> import Data.Map (findWithDefault,fromListWith,toList)

基因的名称可以是任何字符串，对于给定的 Gene g，类型为 PosReg 的函数应该给我们所有g正向调节的基因。

> type Gene = String
> type PosReg = Gene -> [Gene]

从您问题中指定的图表中，我们需要基因的三元组，以便“正调节”关系是传递的，并且子图描述了所需的属性。首先，从图中选择一个任意基因x。接下来，选择 x 正向调节的基因 y 之一。为了保持传递性，z 必须是 x 和 y 正向调节的基因。

> subgraphs :: String -> [(Gene,Gene,Gene)]
> subgraphs g = do
>   x <- choose
>   y <- posRegBy x
>   z <- posRegBy y
>   guard $ z `elem` posRegBy x
>   return (x,y,z)
>   where (choose,posRegBy) = decode g

通过decode中的简单解析器，我们提取了图中的基因列表和一个PosReg函数，该函数给出了受其他基因正向调节的所有基因。

> decode :: String -> ([Gene], PosReg)
> decode g =
>   let pr = fromListWith (++) $ go (lines g)
>       gs = nub $ concatMap (\(a,b) -> a : b) $ toList pr
>   in (gs, (\x -> findWithDefault [] x pr))
>   where
>     go ((words -> [a, op, b]):ls)
>       | op == "+" = (a,[b]) : go ls
>       | otherwise = go ls
>     go _ = []

最后，主程序将它们粘合在一起。对于找到的每个子图，将其打印到标准输出。

> main :: IO ()
> main = mapM_ (putStrLn . show) $ subgraphs graph
>   where graph = "a + b\n\
>                 \b + c\n\
>                 \c - f\n\
>                 \b - d\n\
>                 \a + c\n\
>                 \f + g\n\
>                 \g + h\n\
>                 \f + h\n"

输出：

("a","b","c")
("f","g","h")