dense subgraph 찾기 구현 (MCODE)

 

   지난 번 글에서 언급했던 알고리즘을 C++로 직접 구현해 보자. 실제로 사용한 예는 다음과 같다.

원래의 network 은 다음과 같다. 우리는 다음과 같은 network 에서 dense 한 sub-graph 즉, edge가 많이 연결된 sub-graph를 뽑아 내어야 한다.

원래의 network.

source network 의 sif 파일은 다음과 같다.

edge_file.sif

위에서 dense한 subgraph 를 지금 설명할 코드로 뽑아 내면 다음과 같다. 

아래의 코드로 찾아낸 cluster1 번.

아래의 코드로 찾아낸 cluster2 번.

아래의 코드로 찾아낸 cluster36 번.

node 를 score로 정렬한 이후 seed로 사용되는데, 높은 점수를 갖는 노드부터 sub-graph (난 계속 cluster 라는 이름을 사용했다. 이 글 이후에도 sub-graph 나 cluster를 사용한다)를 찾기 때문에 cluster 번호가 커질수록 cluster의 density 가 낮아 지는 양상을 보여야 한다. 위의 그림에서 보면 이런 패턴이 보이는 것을 알 수 있다. 위의 그래프는 Cytoscape 으로 그린 것이다. MCODE 를 이용하여 위에서 source 로 사용한 네트워크에서 dense 한 cluster를 찾아 보면 default parameter를 사용할 경우 7 개가 찾아 진다. 이것은 각 sub-graph 간에 node 가 절대 겹치지 않게 구하기 때문이다. 즉, MCODE는 disjoint dense sub-graph 를 찾게 된다. 나는 바로 이 disjoint 라는 조건을 버리고, 각 cluster 간에 노드가 중복되도록 해야 했기 때문에 결국 직접 구현하게 되었다.

  전체적인 구조는 다음과 같다.

 

find_complexes 함수가 main 함수이다. 일단 각 node를 scoring 을 한 후, 이 점수를 기준으로 dense 한 지역을 찾게 된다. scoring 을 할 때 k-core 를 이용하는데, 우선 k-core 가 되는 k에 대하여 가장 큰 k 를 찾아야 한다. 이 때, 원 논문에서는 별 말이 없는데, 우리는 k 의 최대값이 node 의 수가 되고, 최소값이 3이라는 것을 알고 있으므로 binary search 를 할 수 있게 된다. 따라서 아래의 코드에서 find_highest_core_k_density 함수에서는 binary search 를 하고 있다. 이제 각 함수를 상위 레벨에서 살펴 보자.

 

find_complexes

// 모든 node에 대하여 점수를 계산한다.

// 노드를 점수에 따라서 내림차순으로 정렬한다.

// 각 노드를 seed로 해서 dense_subgraph 를 찾는다 –> find_complex

find_complex

// 현재 노드가 이미 처리되었으면 그냥 나간다.

// 현재 노드를 처리된 노드 목록에 추가한다.

// 현재 노드의 점수를 얻어 온다. --> score_node

// 현재 노드의 neighborhood nodes NB  얻어 온다.

// NB에 속한 각 노드 n 에 대해서, 현재 노드의 점수 s 와 n의 점수 ns 가 주어진 조건을 만족시키면 n 에 대하여 find_complex를 다시 호출한다.

score_node

// neighbor 가 없는 노드는 점수가 0 이다.

// 현재 노드가 포함된 k-core 중 최대 k 값과 그 k-core의 density d 를 얻어 온다.--> find_highest_core_k_density

// 점수를 k*d 로 한다.

find_highest_core_k_density

// 최소 k인 k_min 을 2로, 최대 k인 k_max 를 0 으로 설정. k_max 의 값을 찾을 max_k 를 0로 설정.  

// edge 개수가 가장 큰 node 의 edge 수를 k_max 로 설정한다.

// 만약 k_max <= k_min 이면 k 를 0, density 를 0.0 으로 하고 종료한다. 

// k_min 과 k_max 가 다른동안 다음을 반복한다.

// -- k = < (k_min + k_max) / 2 > 로 한다. ( 이 때 < a > 는 a를 반올림한 값이다)

// -- 그래프에서 k-core 를 찾는다. --> find_k_core_density

// -- -- 만약 찾아지면 : k_min 을 k 로 설정하고, highest k-core 의 k로 반환시킬 max_k 를 k로 설정, 반환시킬 density  d 를 현재 계산된 density 로 설정.

// -- -- 만약 찾아지지 않으면 : k 와 k_max 가 같으면 중단. 아닐 경우 k_max 를 k 로 설정.

// 위의 반복문이 끝나면 max_k 가 the highest k-core의 k, 이 k-core의 density 는 d가 되므로 이 값을 반환.

find_k_core_density

//// 이 함수는 k 와 density를 동시에 구한다.  

// adjacent matrix 를 만든다.

// 각 노드가 몇 개의 edge를 갖고 있는지를 기록한다.

// 제거된 노드 수 removed_node_num 을 0 으로 설정한다.

// 제거된 노드 수가 전체 노드 수랑 같지 않을동안 다음을 반복.

// -- 모든 노드에 대해서, edge 수가 1 이상, k 미만이면 이 노드를 삭제. 이 때 삭제는 다음을 의미.

// -- -- removed_node_num 을 1 증가시킴.

// -- -- 현재 노드의 edge 수를 0 으로 설정 (just for clear)

// -- -- 현재 노드와 연결된 모든 노드에 대해서, adjacent matrix 를 재설정.

// -- -- -- 연결된 노드의 edge 수를 1개 줄임. 만약 줄어든 값이 0 이면 이 노드 역시 삭제되는 것이므로 removed_node_num 을 1 증가시킴.

// density 계산.

// 삭제된 노드 수가 전체 노드 수보다 작으면 true, 아니면 false 반환.

이제 실제 코드를 보자.

find_complexes

bool find_complexes(bool **graph, int node_number, double d, std::vector<std::set<int> >* dest)

{

    double *score = NULL;

    score_node(graph, node_number, &score);

 

    std::vector<std::pair<int,double> > node2score;

    int i = 0;

    for(i = 0; i<node_number; i++){

        node2score.push_back(std::make_pair<int,double>(i,score[i]));

    }

    std::sort(node2score.begin(), node2score.end(), compare_by_second);

 

    std::set<int> processed_node;

    std::vector<std::pair<int,double> >::const_iterator pos = node2score.begin();

    for(; pos != node2score.end(); pos++){

        const int& node = pos->first;

        if(processed_node.count(node) == 0){

            std::set<int> component;

            find_complex(graph, node_number,score, d, processed_node, node, &component);

 

            if(component.size() > 10){

                dest->push_back(component);

            }

        }

    }

 

    delete score;

 

    return true;

}

위 코드는 역시나 직관적이다^[각주:1]. 위의 pseudo code 와 비교해 가면서 코드를 보면 별로 무리없이 이해할 수 있다.

find_complex

bool find_complex(bool **graph, int node_number, double* score, double d, std::set<int>& processed_node, const int& node, std::set<int> *dest)

{

    if(processed_node.count(node) == 1) return true;

    processed_node.insert(node);

 

    double node_score = score[node];

    std::vector<int> neighbor;

    if(get_neighbor(node, graph, node_number, &neighbor) == false) return false;

    dest->insert(node);

    std::vector<int>::const_iterator n = neighbor.begin();

    for(; n != neighbor.end(); n++){

        double n_score = score[*n];

        if(n_score > (node_score*(1.0 - d))){

            dest->insert(*n);

            find_complex(graph, node_number, score, d, processed_node, *n, dest);

        }

    }

    return true;

}

score_node

bool score_node(bool **graph, int node_number, double **score)

{

    *score = new double[node_number];

    int i = 0;

    for(i = 0; i<node_number; i++){

        (*score)[i] = 0;

        std::vector<int> neighbor;

        if(get_neighbor(i, graph, node_number, &neighbor) == false) continue;

        int k = 0; double density = 0;

        if(find_highest_core_k_density(graph,node_number,neighbor, &k, &density) == false) continue;

        (*score)[i] = k*density;

    }

 

    return true;

}

실제로 부여되는 score는 the highest k-core 에서의 k 와 그 k-core의 density 인 d 를 곱한 값으로 부여하고 있음을 알 수 있다.

find_highest_core_k_density

bool find_highest_core_k_density(bool **graph, int node_number, const std::vector<int>& subgraph_node, int *_k, double *_density)

{

    double density = 0;

    double d = 0;

    int max_k = 0;

    int k_min = 2;

    int k_max = 0;

 

    int all_node_number = (int)(subgraph_node.size());

    // The largest k of k-core in subgraph is upper bounded by

    // the maximum of edge number to a node.

    int i = 0; int j = 0;

    for(i = 0; i<all_node_number; i++){

        int edge_num = 0;

        for(j = 0; j<all_node_number; j++){

            edge_num += (graph[subgraph_node[i]][subgraph_node[j]] == true ? 1 : 0);

        }

        if(k_max < edge_num){

            k_max = edge_num;

        }

    }

 

    if(k_max <= k_min){

        *_k = 0;

        *_density = 0;

        return true;

    }

 

    while(k_min != k_max){

        int k = (int)((k_min + k_max ) * 0.5 + 0.5);

        if(find_k_core_density(graph, subgraph_node, k, &d) == true){

            k_min = k;

            max_k = k;

            density = d;

        }

        else{

            if(k == k_max) break;

            k_max = k;

        }

    }

    *_k = max_k;

    *_density = density;

 

    return true;

}

최대 k-core 를 찾는 것은 원래 O(N³) 이지만 위처럼 하면 그나마 O(N²logN) 으로 줄일 수 있기는 하다.

find_k_core_density

// [graph] is for total graph, and subgraph with nodes in [neighbor] will be

// used to find k-core graph.

bool find_k_core_density(bool **graph, const std::vector<int>& subgraph_node, int k, double* density)

{

    int i = 0; int j = 0;

    int all_node_size = (int)(subgraph_node.size());

    bool **adjacent_matrix = new bool*[all_node_size];

    int *edge_sum = new int[all_node_size];

    for(i = 0; i<all_node_size; i++){

        adjacent_matrix[i] = new bool[all_node_size];

        edge_sum[i] = 0;

        for(j = 0; j<all_node_size; j++){

            adjacent_matrix[i][j] = false;

        }

    }

 

    for(i = 0; i<all_node_size-1; i++){

        for(j = i; j<all_node_size; j++){

            adjacent_matrix[i][j] = adjacent_matrix[j][i] = graph[subgraph_node[i]][subgraph_node[j]];

        }

    }

 

    int removed_node_sum = 0;

    for(i = 0; i<all_node_size; i++){

        int sum = 0;

        for(j = 0; j<all_node_size; j++){

            sum += (adjacent_matrix[i][j] == true ? 1 : 0);

        }

        edge_sum[i] = sum;

        if(sum == 0){ // Is it possible?

            removed_node_sum++;

        }

    }

 

    int removed_node_number = 0;

    do{

        removed_node_number = 0;

        for(i = 0; i<all_node_size; i++){

            if(edge_sum[i] == 0) continue;

            if(edge_sum[i] < k){

                removed_node_number++;

                edge_sum[i] = 0;

                for(j = 0; j<all_node_size; j++){

                    if(i == j) continue;

                    if(adjacent_matrix[i][j] == true){

                        adjacent_matrix[j][i] = false;

                        adjacent_matrix[i][j] = false;

                        edge_sum[j]--;

                        if(edge_sum[j] == 0){

                            removed_node_number++;

                        }

                    }

                }

            }

        }

        removed_node_sum+=removed_node_number;

        if(removed_node_sum > all_node_size){

            std::cout << "ERROR" << std::endl;

            exit(-1);

        }

        if(removed_node_sum == all_node_size) break;

    }

    while(removed_node_number > 0);

 

    int remained_edge_count = 0;

    int remained_node_count = 0;

    for(i = 0; i<all_node_size; i++){

        int edge_count = 0;

        for(j = 0; j<all_node_size; j++){

            if(adjacent_matrix[i][j] == true){

                edge_count++;

            }

        }

        remained_edge_count+=edge_count;

        if(edge_count > 0){

            remained_node_count++;

        }

    }

 

    *density = (double)(remained_edge_count)*2/(remained_node_count*(remained_node_count+1));

 

    for(i = 0; i<all_node_size; i++){

        delete adjacent_matrix[i];

    }

    delete adjacent_matrix;

    delete edge_sum;

 

    return all_node_size - removed_node_sum > 0;

}

이제, 실제로 위의 routine 들을 이용하여 dense한 sub-graph 를 뽑아 내는 루틴을 살펴 보자.

int main(int argc, const char* argv[])

{

    if(argc == 1){

        std::cout << "./dense_graph [edge_file_name]" << std::endl;

        return -1;

    }

 

    std::string edge_file_name = argv[1];

 

 

    bool **graph = NULL;

    std::map<std::string,int> node2index;

    std::map<int,std::string> index2node;

    int node_num = 0;

    int edge_num = 0;

    if(load_graph(edge_file_name, &graph, &node2index, &index2node, &node_num, &edge_num) == false) return -1;

 

    std::cout << edge_file_name << '\t' << node_num << '\t' << edge_num << std::endl;

    int node_number = (int)(node2index.size());

    std::vector<std::set<int> > cluster;

    find_complexes(graph, node_number, 0.05, &cluster); // 0.05가 control 변수.

    if(cluster.empty() == true) return -1;

 

 

    std::ofstream out((edge_file_name + ".cluster.txt").c_str());

    out << "INFO:SOURCE_FILE:\t" << edge_file_name << std::endl;

 

    int i = 1;

    std::vector<std::set<int> >::const_iterator cpos = cluster.begin();

    for(; cpos != cluster.end(); cpos++){

        std::stringstream cluster_name;

        cluster_name << edge_file_name << "_Cluster" << i << ".sif";

        std::ofstream cf(cluster_name.str().c_str());

 

        std::vector<int> c(cpos->begin(), cpos->end());

        std::vector<int>::const_iterator n1 = c.begin();

        int edge_count = 0;

        std::stringstream cout;

        for(; n1 != c.end()-1; n1++){

            std::vector<int>::const_iterator n2 = n1 + 1;

            for(; n2 != c.end(); n2++){

                if(graph[*n1][*n2] == true){

                    edge_count++;

                    cf << index2node[*n1] << "\tpp\t" << index2node[*n2] << std::endl;

                }

            }

        }

        for(n1 = c.begin(); n1 != c.end(); n1++){

            cout << '\t' << index2node[*n1];

        }

 

        int n = (int)(cpos->size());

        double density = (double)(2.0*edge_count)/(n*(n-1));

        out << "Cluster" << i << '\t' << density << '\t' << n << '\t' << edge_count << cout.str() << std::endl;

        i++;

    }

 

    for(i = 0; i<node_number; i++){

        delete graph[i];

    }

    delete graph;

 

    return 0;

}

조금 복잡할 수 있는데, 위 실행 파일이 입력으로 받는 파일은 다음과 같은 구조를 갖고 있다.

EDGE    06436:EDD1    07443:AZIN1    0.996634    0.999722    0.98847    0.999654   
EDGE    00549:EGR1    01270:FOSB    0.999896    0.999962    0.99939    0.999643   
EDGE    01473:PROS1    03443:MEIS2    0.997933    0.99907    0.963672    0.908679
EDGE    00549:EGR1    03881:EGR3    0.996599    0.999486    0.999541    0.998022
...

예제 파일은 다음에 있다.

edge_file.txt

각 line 은 edge 정보를 갖고 있다. 두 번째 column 과 세 번째 column 이 각각의 node 이고, 그 두의 숫자는 어떤 계산된 값^[각주:2]인데 일단 무시하자. 위에서 load_graph 는 위의 edge_file.txt 파일을 읽어서 bool ** 형태의 adjacency matrix 를 채우게 된다. 단, 이 때 위처럼 입력 파일이 text 로 되어 있기 때문에 이것을 숫자로 바꾸는 작업을 하게 된다. 내부적으로 모든 계산은 숫자로 하게 되고, 출력을 할 때만 다시 그 숫자에서 위의 문자열들을 뽑아 내는 형태로 되어 있다. 위의 모든 내용을 담은 cpp 파일은 다음에 있다.

dense_graph2.cpp

컴파일은 간단하다. 필요한 부분을 저 파일에 모두 몰아 넣었기 때문에 간단히 다음과 같이 컴파일 하면 된다.

[adnoctum@bioism dense_graph]$ g++ dense_graph2.cpp -o dense_graph2
[adnoctum@bioism dense_graph]$

dense_graph2.cpp 파일의 95 번째 줄부터가 실제로 위에 필요한 내용들이다^[각주:3]. 위처럼 하고, 다음과 같이 실행시키면 각 cluster 가 sif 파일로 저장이 된고, edge_file.txt.cluster.txt 에는 각 cluster 의 정보들이 들어가 있게 된다. 다음과 같다.

[adnoctum@bioism dense_graph]$ ll
total 380
-rwxrwxr-x 1 adnoctum adnoctum 126101 Dec  3 19:54 dense_graph2
-rw-rw-r-- 1 adnoctum adnoctum  11516 Dec  3 19:19 dense_graph2.cpp
-rw-rw-r-- 1 adnoctum adnoctum  33507 Dec  3 19:33 edge_file.sif
-rw-rw-r-- 1 adnoctum adnoctum 201552 Dec  3 19:10 edge_file.txt
[adnoctum@bioism dense_graph]$ ./dense_graph2 edge_file.txt
edge_file.txt   138     1256
[adnoctum@bioism dense_graph]$ ll
total 560
-rwxrwxr-x 1 adnoctum adnoctum 126101 Dec  3 19:54 dense_graph2
-rw-rw-r-- 1 adnoctum adnoctum  11516 Dec  3 19:19 dense_graph2.cpp
-rw-rw-r-- 1 adnoctum adnoctum  33507 Dec  3 19:33 edge_file.sif
-rw-rw-r-- 1 adnoctum adnoctum 201552 Dec  3 19:10 edge_file.txt
-rw-rw-r-- 1 adnoctum adnoctum   1469 Dec  3 19:57 edge_file.txt_Cluster10.sif
-rw-rw-r-- 1 adnoctum adnoctum   1279 Dec  3 19:57 edge_file.txt_Cluster11.sif
-rw-rw-r-- 1 adnoctum adnoctum   3263 Dec  3 19:57 edge_file.txt_Cluster12.sif
-rw-rw-r-- 1 adnoctum adnoctum   3515 Dec  3 19:57 edge_file.txt_Cluster13.sif
-rw-rw-r-- 1 adnoctum adnoctum   4793 Dec  3 19:57 edge_file.txt_Cluster14.sif
-rw-rw-r-- 1 adnoctum adnoctum   3594 Dec  3 19:57 edge_file.txt_Cluster15.sif
-rw-rw-r-- 1 adnoctum adnoctum   1337 Dec  3 19:57 edge_file.txt_Cluster16.sif
-rw-rw-r-- 1 adnoctum adnoctum   1176 Dec  3 19:57 edge_file.txt_Cluster17.sif
-rw-rw-r-- 1 adnoctum adnoctum   5156 Dec  3 19:57 edge_file.txt_Cluster18.sif
-rw-rw-r-- 1 adnoctum adnoctum    847 Dec  3 19:57 edge_file.txt_Cluster19.sif
-rw-rw-r-- 1 adnoctum adnoctum   4429 Dec  3 19:57 edge_file.txt_Cluster1.sif
-rw-rw-r-- 1 adnoctum adnoctum   3536 Dec  3 19:57 edge_file.txt_Cluster20.sif
-rw-rw-r-- 1 adnoctum adnoctum   7499 Dec  3 19:57 edge_file.txt_Cluster21.sif
-rw-rw-r-- 1 adnoctum adnoctum   1482 Dec  3 19:57 edge_file.txt_Cluster22.sif
-rw-rw-r-- 1 adnoctum adnoctum   4018 Dec  3 19:57 edge_file.txt_Cluster23.sif
-rw-rw-r-- 1 adnoctum adnoctum   1201 Dec  3 19:57 edge_file.txt_Cluster24.sif
-rw-rw-r-- 1 adnoctum adnoctum   1865 Dec  3 19:57 edge_file.txt_Cluster25.sif
-rw-rw-r-- 1 adnoctum adnoctum   1259 Dec  3 19:57 edge_file.txt_Cluster26.sif
-rw-rw-r-- 1 adnoctum adnoctum   1588 Dec  3 19:57 edge_file.txt_Cluster27.sif
-rw-rw-r-- 1 adnoctum adnoctum   1469 Dec  3 19:57 edge_file.txt_Cluster28.sif
-rw-rw-r-- 1 adnoctum adnoctum    594 Dec  3 19:57 edge_file.txt_Cluster29.sif
-rw-rw-r-- 1 adnoctum adnoctum   2103 Dec  3 19:57 edge_file.txt_Cluster2.sif
-rw-rw-r-- 1 adnoctum adnoctum    894 Dec  3 19:57 edge_file.txt_Cluster30.sif
-rw-rw-r-- 1 adnoctum adnoctum   2833 Dec  3 19:57 edge_file.txt_Cluster31.sif
-rw-rw-r-- 1 adnoctum adnoctum    790 Dec  3 19:57 edge_file.txt_Cluster32.sif
-rw-rw-r-- 1 adnoctum adnoctum    864 Dec  3 19:57 edge_file.txt_Cluster33.sif
-rw-rw-r-- 1 adnoctum adnoctum   1137 Dec  3 19:57 edge_file.txt_Cluster34.sif
-rw-rw-r-- 1 adnoctum adnoctum    664 Dec  3 19:57 edge_file.txt_Cluster35.sif
-rw-rw-r-- 1 adnoctum adnoctum    536 Dec  3 19:57 edge_file.txt_Cluster36.sif
-rw-rw-r-- 1 adnoctum adnoctum   2323 Dec  3 19:57 edge_file.txt_Cluster3.sif
-rw-rw-r-- 1 adnoctum adnoctum   3126 Dec  3 19:57 edge_file.txt_Cluster4.sif
-rw-rw-r-- 1 adnoctum adnoctum   6002 Dec  3 19:57 edge_file.txt_Cluster5.sif
-rw-rw-r-- 1 adnoctum adnoctum   3263 Dec  3 19:57 edge_file.txt_Cluster6.sif
-rw-rw-r-- 1 adnoctum adnoctum   1196 Dec  3 19:57 edge_file.txt_Cluster7.sif
-rw-rw-r-- 1 adnoctum adnoctum   1518 Dec  3 19:57 edge_file.txt_Cluster8.sif
-rw-rw-r-- 1 adnoctum adnoctum   7844 Dec  3 19:57 edge_file.txt_Cluster9.sif
-rw-rw-r-- 1 adnoctum adnoctum   8773 Dec  3 19:57 edge_file.txt.cluster.txt
[adnoctum@bioism dense_graph]$

standard out 으로 출력된 위의 붉은 글씨는, 입력 파일과 node 수, edge의 수이다. 그냥 보려고 출력했다.



아... 써 놓고 보니 조금 복잡해 보인다. 아마도 pseudo code 과 직접적으로 비교해 보면 실제 코드의 이해가 어느 정도 수월하리라 생각한다. 나는 될 수 있으면 직관적으로 코딩하려 하기 때문에... 만약 위의 실제 코드 내용 중 난해한 부분에 대한 지적이 있으면 좀 더 자세한 설명을 추가한다.

1. 나는 어지간히 시간을 잡아 먹는 코드가 아닌 경우에는 위처럼 매우 직관적으로 코드를 짠다. 다소 비효율적으로 보이는 것 같은 것을 알고 있는 경우라 하더라도 직관적인 코드가 가능하면 그렇게 한다. [본문으로]
2. 실제로는 pearson's correlation coefficient 에 대한 rank 이다. 이것은 내 연구주제에 관련된 것이기 때문에 여기서 말할 필요는 없으므로 그냥 넘어 간다. [본문으로]
3. dense_graph2.cpp 에서 왜 2 라는 숫자가 붙었느냐 하면, 원래 나는 그래프를 multimap 으로 표현하고 그것이 dense_graph 였는데, 너무 느려서 위처럼 이중포인터를 이용한 자료형으로 바꾸면서 2 를 붙이게 되었다. [본문으로]