こんにちは、受託コンサルティングチームの hosor です。
私は技術ブログでアメリエフを知ったので、執筆する側になる日が来て感慨深いです。
以前に、以下のようなブログを投稿していました。
こちらでは、着目パスウェイの第一階層に該当するパスウェイを検索していました。
Reactome Pathwayの最上位カテゴリを自動的に検索する方法 - アメリエフの技術ブログ
今回は、着目パスウェイの第二階層に該当するパスウェイの検索方法について書いていこうと思います。
得られたパスウェイを第二階層でまとめることで、全体の傾向を把握しつつ biological な観点からの解釈も可能になります。
作業環境
- CentOS 7
- Python v3.11.3
- Python 標準ライブラリ collections
使用データ
Reactome の Download より、以下の 2 つのデータを取得しました。
今回は、前者を pathways_hierarchy_relationship.txt と、後者を complete_list_of_pathways.txt と命名しています。
スクリプト
まず、着目パスウェイの 1 つ上の階層のパスウェイを検索できるようにします。
同時に、すべてのパスウェイを all_pathways に入れていきます。
from collections import defaultdict
upper_level = defaultdict(list)
all_pathways = set()
pathways_hierarchy_relationship_input = open("pathways_hierarchy_relationship.txt", 'r')
for row in pathways_hierarchy_relationship_input:
lst = list(map(str, row.split()))
# パスウェイ lst[1] の 1 つ上層のパスウェイは lst[0]
upper_level[lst[1]].append(lst[0])
all_pathways.add(lst[0])
all_pathways.add(lst[1])
次に、着目パスウェイを第二階層のパスウェイに変換できるようにします。
ここでは、着目パスウェイから、第二階層に辿り着くまで上層のパスウェイに上がっていく、ということを行なっています。
注意点
着目パスウェイが第一階層のパスウェイである場合は変換ができないため Unidentified と変換するようにしています。
着目パスウェイに第二階層が複数存在するケース (R-HSA-354194) や、重複するケース (R-HSA-166058) があるため、set で管理しています。
from collections import deque
target2second_level = defaultdict(set)
for target_pathway in all_pathways:
# 注目しているパスウェイが第一階層である場合
if not target_pathway in upper_level:
target2second_level[target_pathway].add("Unidentified")
continue
que = deque()
que.append(target_pathway)
while len(que) > 0:
cur_pathway = que.popleft()
for cur_upper_level in upper_level[cur_pathway]:
if cur_upper_level in upper_level:
# 1 つ上の階層が第一階層ではない場合
# 1 つ上の階層に移動する
que.append(cur_upper_level)
else:
# 1 つ上の階層が第一階層である場合
# cur_pathway は第二階層である
target2second_level[target_pathway].add(cur_pathway)
以下のように、正しく検索することができました！
print(target2second_level["R-HSA-354194"])
# {'R-HSA-354192', 'R-HSA-76002'}
print(target2second_level["R-HSA-166058"])
# {'R-HSA-168249'}
また、ID を Description に変換することも可能です。
ID2Description = defaultdict(str)
complete_list_of_pathways_input = open("complete_list_of_pathways.txt", 'r')
for row in complete_list_of_pathways_input:
lst = list(map(str, row.split('\t')))
ID2Description[lst[0]] = lst[1]
print([ID2Description[list(target2second_level["R-HSA-354194"])[0]], ID2Description[list(target2second_level["R-HSA-354194"])[1]]])
# ['Integrin signaling', 'Platelet activation, signaling and aggregation']
print([ID2Description[list(target2second_level["R-HSA-166058"])[0]]])
# ['Innate Immune System']
まとめ
Reactome のデータを使うことで、第二階層を検索することができる！
