前回の記事でGO解析ツール「goatools」について紹介しました。
その時のコマンドライン操作のみからのGO関係可視化図の見栄えがいまいちだったので、見栄えがよくなるように改良を試みました。
改良による見栄えの改善に関する備忘録を残しておきます。

実行環境設定

goatoolsのインストールについては前回の記事を参照

色情報を取り扱うライブラリcolourが前回から追加で必要
pip install colour

Ubuntu20.04 & Python3.8.10 にて実行環境構築を実施

前回までのGO関係可視化の問題点

goatoolsのコマンドツール「go_plot.py」により下記のGO関係可視化図が前回得られた。

echo -e "#ff2200\tGO:0009409\n#ff4400\tGO:0009628\n#ff6600\tGO:0006950\n#ff8800\tGO:0009266\n#ffaa00\tGO:0009415\n#ffbb00\tGO:0001101\n#ffcc00\tGO:0050896\n#ffdd00\tGO:0009414\n#ffee00\tGO:0010035\n#ffff00\tGO:0006970" \
> go_enrichment_color_list.txt
go_plot.py --go_file=go_enrichment_color_list.txt --outfile=plot_color_go_enrichment_BP.png

GO関係可視化図を見ると見栄えに関して以下の点が気になり、何とかしたいなと感じました。

GO機能説明が横に長いため、図全体が横長で見づらい
「L1 D1 d1979」のようなGOの階層・深さ等の記載が邪魔に感じる
pvalueを元に色の強弱をつける想定で可視化したが、pvalueの記載が図中にないので分かりづらい。

GO関係可視化の改良結果

下記の改良をしたPythonスクリプトを作成しました(スクリプトは一番最後に記載)。

GO機能説明が20文字以上の場合に改行
GOの階層・深さ等の出力をしないように設定
pvalueをGO機能説明とともに出力する
pvalue値の大小関係を表現する段階色を自動選択する

見栄えを改良したコード実行

python color_go_plot.py color_go_plot.png

出力結果図を見ると、問題点が改良され、見栄えが良くなった。

その他例：可視化するGO数を増やして段階色の指定を変更
(ソース内での段階色指定について、convert_hexcolor_gradient()関数の引数にてfrom_color="lightblue", to_color="green"に変更)

実行ソース

とりあえずのサンプルコードなので、入力するGOとPvalueは現状ではコード内に埋め込んでいる。
外部ファイルからGOと対応するPvalueを入力できるようにすれば、より汎用的に利用できるようになる。
このソースをベースにし、目的に応じてコードを書き換えれば、ある程度のプロット作成に対応できそう。

color_go_plot.py

import math
import sys
from typing import Dict, List, Optional

import numpy as np
from colour import Color
from goatools.godag_obosm import OboToGoDagSmall
from goatools.godag_plot import GODagPltVars, GODagSmallPlot
from goatools.obo_parser import GODag


def main(plot_outfile: str):
    """Plot Color GO DAG(Directed Acyclic Graph)

    Args:
        plot_outfile (str): Plot output file path
    """
    # GOenrichment result example
    goid2pvalue = {
        "GO:0009409": 4.42507619016377e-16,
        "GO:0009628": 1.70399392910078e-14,
        "GO:0006950": 1.70399392910078e-14,
        "GO:0009266": 2.58192361031249e-14,
        "GO:0009415": 6.399674360673e-14,
        "GO:0001101": 6.44890640551705e-14,
        "GO:0050896": 3.4574338076728e-13,
        "GO:0009414": 3.81276641723126e-13,
        "GO:0010035": 4.04311940400782e-11,
        "GO:0006970": 1.48573634972923e-07,
        # "GO:0009737": 9.21442857296839e-07,
        # "GO:0097305": 9.35271973837119e-07,
        # "GO:0042221": 9.35271973837119e-07,
        # "GO:1901700": 3.44241932707208e-06,
        # "GO:0009269": 7.13392304953131e-06,
        # "GO:0009631": 1.19992021145805e-05,
        # "GO:0033993": 2.06259190472414e-05,
        # "GO:0009725": 0.000458054576231,
        # "GO:0009719": 0.000458054576231,
    }

    # Convert pvalue to common logarithm(log10) of the absolute value
    pvalue_abs_log10_list = [abs(np.log10(v)) for v in goid2pvalue.values()]

    # Get hexcolor from converted pvalue
    pvalue_hexcolor_list = convert_hexcolor_gradient(
        pvalue_abs_log10_list, from_color="yellow", to_color="red"
    )
    goid2color = {
        go: hexcolor for go, hexcolor in zip(goid2pvalue.keys(), pvalue_hexcolor_list)
    }

    # Plot color go dag
    color_go_plot(plot_outfile, goid2color, goid2pvalue)


def convert_hexcolor_gradient(
    value_list: List[float],
    step_num: int = 20,
    default_min: Optional[float] = 0,
    default_max: Optional[float] = 10,
    from_color: str = "yellow",
    to_color: str = "red",
) -> List[str]:
    """Convert to gradient hexcolor (e.g. "#ffff00") list based on the size of the value

    Args:
        value_list (List[float]): List of float values
        step_num (int, optional): Number of gradient steps. Defaults to 10.
        default_min (float, optional): Default minimum value. Defaults to 0.05.
        default_max (float, optional): Default maximum value. Defaults to 10.
        from_color (str, optional): Start color for the gradient. Defaults to "yellow".
        to_color (str, optional): End color of the gradient. Defaults to "red".

    Returns:
        List[str]: List of gradient hexcolor
    """
    # Create gradient hexcolor list
    color_gradient_list = list(Color(from_color).range_to(Color(to_color), step_num))
    hexcolor_gradient_list = [color.get_hex_l() for color in color_gradient_list]

    # Define value range for gradient hexcolor
    min_value, max_value = min(value_list), max(value_list)
    if default_min:
        min_value = min(min_value, default_min)
    if default_max:
        max_value = max(max_value, default_max)
    value_range_list = list(np.linspace(min_value, max_value, step_num + 1))[1:]

    # Get gradient hexcolor corresponding to input value
    convert_hexcolor_gradient_list = []
    for value in value_list:
        for idx, value_range in enumerate(value_range_list):
            if value <= value_range:
                convert_hexcolor_gradient_list.append(hexcolor_gradient_list[idx])
                break

    return convert_hexcolor_gradient_list


def color_go_plot(
    plot_outfile: str,
    goid2color: Dict[str, str],
    goid2pvalue: Dict[str, float] = {},
    obo_file: str = "go-basic.obo",
) -> None:
    """Plot GO DAG using self-defined GO color

    Args:
        plot_outfile (str): Output plot file path
        goid2color (Dict[str, str]): go id and hexcolor dict
        goid2pvalue (Dict[str, float], optional): go id and pvalue dict. Defaults to {}.
        obo_file (str, optional): OBO file path. Defaults to "go-basic.obo".
    """
    # Get plot target GO DAG
    obodag = GODag(obo_file)
    godagsmall = OboToGoDagSmall(goids=goid2color.keys(), obodag=obodag).godag

    # Wrapping GO description line at appropriate location
    for v in godagsmall.go2obj.values():
        if len(v.name) < 20:
            continue
        split_word = v.name.split(" ")
        split_cnt = math.ceil(len(split_word) / 2) - 1
        line_wrap_name = ""
        for cnt, word in enumerate(split_word):
            newline_or_space = "\n" if cnt == split_cnt else " "
            line_wrap_name += word + newline_or_space
        v.name = line_wrap_name

    # Add pvalue to the end of GO description
    for v in godagsmall.go2obj.values():
        if v.id in goid2pvalue.keys():
            v.name += f"\n{goid2pvalue[v.id]:.2e}"

    # Suppress useless header plot (e.g. L2 D2)
    godag_plg_vars = GODagPltVars()
    godag_plg_vars.fmthdr = "{GO}"

    # Create plot obj & add plot color
    godagplot = GODagSmallPlot(godagsmall, abodag=obodag, GODagPltVars=godag_plg_vars)
    godagplot.goid2color = goid2color

    # Plot color go dag
    godagplot.plt(plot_outfile, "pydot")


if __name__ == "__main__":
    args = sys.argv
    main(args[1])

参考

goatoolsによるGO解析方法まとめ