5月 6, 2022 / 最終更新日時 : 5月 6, 2022 Shuji Suzuki (shu) プログラミング

ReactomeからPathwayの階層構造とPathwayに関連するGeneのデータを取得する

1 はじめに
2 P-Netで利用されているReactomeのデータ
3 Reactomeから大本のデータを取得する
4 Reactomeから落としたデータをPythonで読み込み
5 Reactomeのデータを組み合わせて使う
6 終わりに

はじめに

だいぶ前ですが以下の論文を読み、ReactomeからどのようにしてPathwayデータを取得するのか？というのが気になっていました。

Elmarakeby, H.A., Hwang, J., Arafeh, R. et al. Biologically informed deep neural network for prostate cancer discovery. Nature 598, 348–352 (2021). https://doi.org/10.1038/s41586-021-03922-4

こちらの論文自体はReactomeというPathwayデータベースのデータを組み合わせてP-NETという深層学習モデルを構築して前立腺がんの患者のデータに対して適用した研究になります。

この深層学習のモデルを構築する際、ReactomeからPathwayに関するデータを取得して利用しているのですが、ReactomeからどうPathwayのデータを取得して、それをどう加工すればよさそうか？が論文を読んだだけでは分からず、コードを読んで調べたので本日の記事はそのまとめになります。

今回の記事で利用したコードはこちらに置いてあります。

https://github.com/shu65/reactome-example/blob/main/Reactome_gene_pathway_hierarchy_relationship.ipynb

P-Netで利用されているReactomeのデータ

P-NetではPathwayの階層情報に基づいて深層学習のモデルを構築しています。この際、利用されているデータは以下のものになります。

Pathway毎に関連するGeneの集合
Pathwayの階層構造
Pathwayの名前と種

これらをReactomeから取得する方法を紹介していきます。わかりやすいようにReactomeのリンクやページのスクリーンショットも合わせのせています。これらは2022/05/06時点のものになります。アクセスする時期によってはこれらが変わっている可能性もあるので注意してください。

Reactomeから大本のデータを取得する

Reactomeは各種データを様々なフォーマットでダウンロードできるようになっています。それらはここにまとまっています。

https://reactome.org/download-data

今回利用するデータはそれぞれ以下のファイルに書かれた情報から抽出することができます。

Pathway毎に関連するGeneの集合

Pathway毎に関連するGeneの集合のデータは[Specialized data formats]→[Reactome Pathways Gene Set.]から落とすことができます。

こちらはGMTフォーマットで書かれたファイルになります。具体的には行ごとに1つのPathwayについて書かれており、カラムはPathwayの名前、Reactome Stable identifiers (ST_ID)が続き、それ以降は登場するGeneが並んでいます。

Pathwayの階層構造

Pathwayの階層構造のデータは[Pathways] → [Pathways hierarchy relationship]から落とすことができます。

[Read more] のリンクに詳しいフォーマットが書かれていますが、タブ区切りで最初のカラムが親のPathwayのST_ID、2つ目が子のST_IDになります。

Pathwayの名前と種

Pathwayの名前と種に関するデータはPathwayの階層構造のデータは[Pathways] → [Complete List of Pathways]から落とすことができます。

[Read more] のリンクに詳しいフォーマットが書かれていますが、タブ区切りで最初のカラムがPathwayのReactome Stable identifiers (ST_ID)、2つ目がPathwayの名前、3つ目が種になります。

Reactomeから落としたデータをPythonで読み込み

ここまででほしいデータがどこから落とすことができるか説明しました。ここからは実際に使う際に利用しやすいようにPythonで読み込むパーサーを書いたのでその紹介をします。

Pathway毎に関連するGeneの集合のデータのパース

Reactomeから落としてきたGMTファイルをパースするスクリプトは以下の通りです。

import pandas as pd


def read_reactome_gmt(file_path):
  data_dict_list = []
  with open(file_path) as f:
      for i, line in enumerate(f):
          values = line.strip().split("\t")
          st_id = values[1]
          genes = values[3:]
          for gene in genes:
              data_dict_list.append({'st_id': st_id, 'gene': gene})
  df = pd.DataFrame(data_dict_list)
  return df

以下のようにようにファイルのパスを渡すとpandasのデータフレームでファイルから読み込んだ結果を返すようにしてあります。

pathway_gene_df = read_reactome_gmt("ReactomePathways.gmt")
print(pathway_gene_df)

出力は以下の通りです。

               st_id   gene
0       R-HSA-164843  HMGA1
1       R-HSA-164843   LIG4
2       R-HSA-164843  PSIP1
3       R-HSA-164843  XRCC4
4       R-HSA-164843  XRCC5
...              ...    ...
121252  R-HSA-192905     NP
121253  R-HSA-192905     NS
121254  R-HSA-192905     PA
121255  R-HSA-192905    PB1
121256  R-HSA-192905    PB2

[121257 rows x 2 columns]

Pathwayの階層構造の読み込み

Reactomeから落としたPathwayの階層構造のデータを読み込むパーサーは以下の通りです。

import pandas as pd


def read_reactome_pathway_hierarchy_relationship(file_path):
  df = pd.read_csv(file_path, sep='\t')
  df.columns = ['parent_st_id', 'child_st_id']
  return df

Reactomeから落としたファイルパスを渡すと以下のようなpandasのデータフレームを返します。

pathway_hierarchy_df = read_reactome_pathway_hierarchy_relationship("ReactomePathwaysRelation.txt")
print(pathway_hierarchy_df)

出力は以下の通りです。

       parent_st_id    child_st_id
0      R-BTA-109581  R-BTA-5357769
1      R-BTA-109581    R-BTA-75153
2      R-BTA-109582   R-BTA-140877
3      R-BTA-109582   R-BTA-202733
4      R-BTA-109582   R-BTA-418346
...             ...            ...
21521  R-XTR-983705  R-XTR-5690714
21522  R-XTR-983705   R-XTR-983695
21523  R-XTR-983712  R-XTR-2672351
21524  R-XTR-983712   R-XTR-936837
21525  R-XTR-991365   R-XTR-997272

[21526 rows x 2 columns]

Pathwayの名前と種の読み込み

ReactomeのPathwayの名前と種のパーサーは以下の通りです。

import pandas as pd


def read_reactome_complete_list(file_path):
  df = pd.read_csv(file_path, sep='\t')
  df.columns = ['st_id', 'pathway_name', 'species']
  return df

こちらも以下のようにしてファイルパスを指定するとpandasのデータフレームを返します。

complete_list = read_reactome_complete_list("ReactomePathways.txt")
print(complete_list)

complete_list = read_reactome_complete_list("ReactomePathways.txt")
print(complete_list)
               st_id                                       pathway_name  \
0      R-BTA-1971475  A tetrasaccharide linker sequence is required ...   
1      R-BTA-1369062              ABC transporters in lipid homeostasis   
2       R-BTA-382556             ABC-family proteins mediated transport   
3      R-BTA-9033807                       ABO blood group biosynthesis   
4       R-BTA-418592          ADP signalling through P2Y purinoceptor 1   
...              ...                                                ...   
21417   R-XTR-193639                           p75NTR signals via NF-kB   
21418   R-XTR-111995                               phospho-PLA2 pathway   
21419   R-XTR-191859                                     snRNP Assembly   
21420   R-XTR-379724                                tRNA Aminoacylation   
21421   R-XTR-199992                trans-Golgi Network Vesicle Budding   

                  species  
0              Bos taurus  
1              Bos taurus  
2              Bos taurus  
3              Bos taurus  
4              Bos taurus  
...                   ...  
21417  Xenopus tropicalis  
21418  Xenopus tropicalis  
21419  Xenopus tropicalis  
21420  Xenopus tropicalis  
21421  Xenopus tropicalis  

[21422 rows x 3 columns]

Reactomeのデータを組み合わせて使う

これらのデータを組み合わせて利用する例として、ヒトのPathwayの中で一番上の親のPathwayとその一つ下のPathway、2つの階層のPathwayの子のリストを取得し、その子の一つのPathwayの名前とGeneのリストを取得するコードを示します。

まずはヒトのPathwayを抽出してPythonのNetworkXを利用して有向グラフを作り、一番上の親をrootというノードにつなげたグラフを作ります。

import networkx as nx

human_pathway_ids = complete_list[complete_list["species"] == 'Homo sapiens']["st_id"]
human_pathway_hierarchy_df = pathway_hierarchy_df[pathway_hierarchy_df["parent_st_id"].isin(human_pathway_ids) & pathway_hierarchy_df["child_st_id"].isin(human_pathway_ids)]
human_pathway_graph = nx.from_pandas_edgelist(human_pathway_hierarchy_df, source="parent_st_id", target="child_st_id", create_using=nx.DiGraph())
root_pathways = [n for n, d in human_pathway_graph.in_degree() if d==0] 
root_edges = [("root", n) for n in root_pathways] 
human_pathway_graph.add_edges_from(root_edges)

これでヒトのPathwayの階層構造を示した有向グラフができました。これを使えばrootからの最小距離が2以下のnodeを列挙することで、一番上の親のPathwayとその一つ下のPathway、2階層のPathwayが取得できます。

selected_ids = nx.single_source_shortest_path_length(human_pathway_graph, source="root", cutoff=2)

selected_idsをprintするとこのような形になります。

length 1 {'R-HSA-109581': 2,
 'R-HSA-109582': 1,
 'R-HSA-112307': 2,
 'R-HSA-112315': 2,
 'R-HSA-112316': 1,
 'R-HSA-1181150': 2,
 'R-HSA-1187000': 2,
 'R-HSA-1266738': 1,
...

ここからsuccessors()を使って各Pathwayの子のPathwayを取得します。

for pathway_id, shortest_path_length in selected_ids.items():
  if shortest_path_length > 0:
    children = list(human_pathway_graph.successors(pathway_id))
    print("length", shortest_path_length, "parent", pathway_id, "children", children)

ちなみに、shortest_path_length=0はrootだけなのでskipしています。出力としては以下のようになります。

length 1 parent R-HSA-1852241 children ['R-HSA-1592230', 'R-HSA-5617833']
length 1 parent R-HSA-5357801 children ['R-HSA-109581', 'R-HSA-5218859']
length 1 parent R-HSA-1266738 children ['R-HSA-1181150', 'R-HSA-186712', 'R-HSA-381340', 'R-HSA-452723', 'R-HSA-525793', 'R-HSA-5619507', 'R-HSA-5682910', 'R-HSA-6805567', 'R-HSA-9616222', 'R-HSA-9675108', 'R-HSA-9690406']
length 1 parent R-HSA-4839726 children ['R-HSA-3247509']
length 1 parent R-HSA-9709957 children ['R-HSA-2187338', 'R-HSA-381753', 'R-HSA-9659379', 'R-HSA-9717189']
length 1 parent R-HSA-1474244 children ['R-HSA-1474228', 'R-HSA-1474290', 'R-HSA-1566948', 'R-HSA-1566977', 'R-HSA-216083', 'R-HSA-3000157', 'R-HSA-3000171', 'R-HSA-3000178', 'R-HSA-8941237']
length 1 parent R-HSA-9612973 children ['R-HSA-1632852', 'R-HSA-9613829', 'R-HSA-9615710']
length 1 parent R-HSA-397014 children ['R-HSA-390522', 'R-HSA-445355', 'R-HSA-5576891']
...

ここで試しにR-HSA-1592230の名前を出してみます。コードとしては以下のようになります。

print(complete_list[complete_list["st_id"] == "R-HSA-1592230"])

出力は以下の通りです。

               st_id              pathway_name       species
11395  R-HSA-1592230  Mitochondrial biogenesis  Homo sapiens

また、R-HSA-1592230の関連するgeneは以下のようにして出力できます。

print(pathway_gene_df[pathway_gene_df["st_id"] == "R-HSA-1592230"])

出力は以下の通りです。

               st_id    gene
65977  R-HSA-1592230   ALAS1
65978  R-HSA-1592230    APOO
65979  R-HSA-1592230   APOOL
65980  R-HSA-1592230    ATF2
65981  R-HSA-1592230   ATP5B
...              ...     ...
66066  R-HSA-1592230    TFAM
66067  R-HSA-1592230   TFB1M
66068  R-HSA-1592230   TFB2M
66069  R-HSA-1592230    TGS1
66070  R-HSA-1592230  TMEM11

[94 rows x 2 columns]

終わりに

Reactomeで公開されているデータを読み込んでPathwayの階層構造やPathwayに関連するGeneのデータを読み込む方法を紹介しました。最初はGraph Databaseから頑張って抽出しないといけないのかと思っていましたが、調べてみると簡単なことがわかりました。ただ、何も知らない状態ではどうしていいのかわからないことが多かったので記事にまとめてみました。同じように悩んでいる方の参考になれば幸いです。

カテゴリー: プログラミング

タグ: Bioinformatics

ReactomeからPathwayの階層構造とPathwayに関連するGeneのデータを取得する

はじめに

P-Netで利用されているReactomeのデータ

Reactomeから大本のデータを取得する

Pathway毎に関連するGeneの集合

Pathwayの階層構造

Pathwayの名前と種

Reactomeから落としたデータをPythonで読み込み

Pathway毎に関連するGeneの集合のデータのパース

Pathwayの階層構造の読み込み

Pathwayの名前と種の読み込み

Reactomeのデータを組み合わせて使う

終わりに

コメントを残すコメントをキャンセル

[勉強ノート] 「線形計算の数理」 7.7.1 Arnoldi法

UniProtのWeb APIを使ってほしい遺伝子のタンパク質配列を取ってくる

はじめに

P-Netで利用されているReactomeのデータ

Reactomeから大本のデータを取得する

Pathway毎に関連するGeneの集合

Pathwayの階層構造

Pathwayの名前と種

Reactomeから落としたデータをPythonで読み込み

Pathway毎に関連するGeneの集合のデータのパース

Pathwayの階層構造の読み込み

Pathwayの名前と種の読み込み

Reactomeのデータを組み合わせて使う

終わりに

コメントを残す コメントをキャンセル

コメントを残すコメントをキャンセル