Aurora blog

バイオインフォ・バイオテクノロジーにまつわる情報

Knowledge graph①: Neo4j Graph Academy備忘録

背景

グラフデータベースNeo4jに触る機会がありNeo4j Graph Academyを受講した。忘れそうなポイントをこのページにメモとして残す。

Cypher

MATCH/OPTIONAL MATCH

  • OPTIONAL MATCH: 該当のクエリが存在する場合にのみ結果を返す (存在しない場合はnullを返す)
MATCH (m:Movie) WHERE m.title = "Kiss Me Deadly"
MATCH (m)-[:IN_GENRE]->(g:Genre)<-[:IN_GENRE]-(rec:Movie)
OPTIONAL MATCH (m)<-[:ACTED_IN]-(a:Actor)-[:ACTED_IN]->(rec)
RETURN rec.title, a.name

特定のプロパティのみを持つノードを返す場合

MATCH (m:Movie)
WHERE m.title CONTAINS 'Matrix'
RETURN m { .title, .released } AS movie

MERGE/DELETE

新しいノード・エッジを作成/削除する。以下の性質があるので注意。

  • MERGEは条件に該当するノード・エッジが既に存在する場合は実行されない。
  • MERGEはエッジの方向性を指定しないと「左側のノードから右側のノード」という方向性がデフォルトでアサインされる *1
  • DETACH DELETE: DELETEではエッジを持つノードを削除できない様になっている。こういったノードを削除する場合に使う。
// 以下の様に実行すると全てのノード・エッジが消える (=グラフが完全消去される)
MATCH (n)
DETACH DELETE n

WHERE

以下の演算子が存在

  • =: Equal
  • <>: Not equal
  • >/<: Greater / Less than
  • OR: OR演算子
  • A STARTS WITH B: 文字列AがBで始まっているか
  • A ENDS WITH B: 文字列AがBで終わっているか
  • A CONTAINS B: 文字列AにBが含まれているか

プロパティが存在している場合は以下のように検証可能

MATCH (p:Person)
WHERE p.died IS NOT NULL

SET系

MERGE/MATCHで指定したノード・エッジのプロパティを追加する

  • SET: 常に実行される
  • ON CREATE SET: MERGEによってノード・エッジが作られた時にのみ実行される
  • ON MATCH SET: MERGE実行時に既にノード・エッジが存在する場合にのみ実行される
// Find or create a person with this name
MERGE (p:Person {name: 'McKenna Grace'})

// Only set the `createdAt` property if the node is created during this query
ON CREATE SET p.createdAt = datetime()

// Only set the `updatedAt` property if the node was created previously
ON MATCH SET p.updatedAt = datetime()

// Set the `born` property regardless
SET p.born = 2006

RETURN p

SETではプロパティを追加するだけでなくノードにラベル (Node Type) を追加する場合にも使える

MATCH (p:Person)
WHERE exists ((p)-[:DIRECTED]-())
SET p:Director

プロパティを消す場合にも使える

MATCH (m:Movie)
SET m.genres = null

WITH

WITHを使うことでMATCHの前に特定の値で変数を定義することができる

WITH 'Tom Hanks' AS actorName
MATCH (p:Person)-[:ACTED_IN]->(m:Movie)
WHERE p.name = actorName
RETURN m.title AS movies

WITHを挟むとRETURN文に記載できる変数名がリセットされる。以下のケースのようにRETURNで返したい変数 (e.g., p, m) がある場合はWITHに入れる必要がある。

WITH  'Clint Eastwood' AS a, 'high' AS t
MATCH (p:Person)-[:ACTED_IN]->(m:Movie)
WITH p, m, toLower(m.title) AS movieTitle
WHERE p.name = a
AND movieTitle CONTAINS t
RETURN p.name AS actor, m.title AS movie

WITHを使った事例: 平均スコアの計算

MATCH (:Movie {title: 'Toy Story'})-[:IN_GENRE]->(g:Genre)<-[:IN_GENRE]-(m)
WHERE m.imdbRating IS NOT NULL
WITH g.name AS genre,
count(m) AS moviesInCommon,
sum(m.imdbRating) AS total
RETURN genre, moviesInCommon,
total/moviesInCommon AS score
ORDER By score DESC

UNWIND

プロパティのリストに含まれる項目ごとに操作を実行することができる。モデルをリファクタリングする際に使用する事例が紹介されていた。

MATCH (m:Movie)
UNWIND m.genres AS genre
MERGE (g:Genre {name: genre})
MERGE (m)-[:IN_GENRE]->(g)
SET m.genres = null

CREATE CONSTRAINT

特定のプロパティ (e.g., ID) に重複を許さないなどの制約を設けることができる。検索の速度改善にもつながるらしい。

CREATE CONSTRAINT Genre_name IF NOT EXISTS
FOR (x:Genre)
REQUIRE x.name IS UNIQUE

ORDER BY

結果を何らかのプロパティでソートする場合に使う

  • DESC: 降順
MATCH (p:Person)
WHERE p.born IS NOT NULL
RETURN p.name AS name, p.born AS birthDate
ORDER BY p.born DESC

2つ以上の条件を使う場合

MATCH (m:Movie)<-[ACTED_IN]-(p:Person)
WHERE m.imdbRating IS NOT NULL
RETURN m.title, m.imdbRating, p.name, p.born
ORDER BY m.imdbRating DESC, p.born DESC

CASE-WHEN

MATCH (m:Movie)<-[:ACTED_IN]-(p:Person)
WHERE p.name = 'Henry Fonda'
RETURN m.title AS movie,
CASE
WHEN m.year < 1940 THEN 'oldies'
WHEN 1940 <= m.year < 1950 THEN 'forties'
WHEN 1950 <= m.year < 1960 THEN 'fifties'
WHEN 1960 <= m.year < 1970 THEN 'sixties'
WHEN 1970 <= m.year < 1980 THEN 'seventies'
WHEN 1980 <= m.year < 1990 THEN 'eighties'
WHEN 1990 <= m.year < 2000 THEN 'nineties'
ELSE  'two-thousands'
END
AS timeFrame

UNWIND

リストを分解、行へと分ける

MATCH (m:Movie)-[:ACTED_IN]-(a:Actor)
WHERE a.name = 'Tom Hanks'
UNWIND m.languages AS lang
RETURN m.title AS movie,
m.languages AS languages,
lang AS language

CALL

サブクエリをつくる。例えば下記のように実行すればMovieのRatingについての検索範囲を絞る (=効率を上げる) ことが可能。

CALL {
   MATCH (m:Movie) WHERE m.year = 2000
   RETURN m ORDER BY m.imdbRating DESC LIMIT 10
}
MATCH  (:User)-[r:RATED]->(m)
RETURN m.title, avg(r.rating)

以下のように途中にサブクエリを挟むことも可能

MATCH (m:Movie)
CALL {
    WITH m
    MATCH (m)<-[r:RATED]-(u:User)
     WHERE r.rating = 5
    RETURN count(u) AS numReviews
}
RETURN m.title, numReviews
ORDER BY numReviews DESC

UNION

2つのクエリの出力を連結する。2つのクエリは同じ項目をリターンする必要がある。

  • UNION ALL: ALLをつけないと2つのクエリで重複する行は除外される (DISTINCTと同等)が、つけると重複も結果として出力される
MATCH (m:Movie) WHERE m.year = 2000
RETURN {type:"movies", theMovies: collect(m.title)} AS data
UNION ALL
MATCH (a:Actor) WHERE a.born.year > 2000
RETURN { type:"actors", theActors: collect(DISTINCT a.name)} AS data

Methods

  • exists: 引数のクエリに該当するノード・エッジが存在するかTRUE/FALSEを返す
  • not exists: 引数のクエリに該当するノード・エッジが存在しない場合にTRUEを返す
  • keys: ラベル・リレーションに存在するプロパティの一覧を返す
  • coalesce: 第一引数がnullでなければ第一引数を、nullであれば第二引数を返す
  • split: string型を特定の文字でsplitしてlist型へ変換する
  • toInteger: int型へ変換
  • type: ラベル・リレーションのタイプを返す
  • size: リストのサイズを返す
  • date/datetime/time: 現時刻をdate/datetime/time型で返す
  • shortestPath: 最短経路を返す
  • trim: stringの先頭・末尾にある空白を取り除く

以下の関数はAggregation (group by)する

  • count: 行数をカウントする
  • collect: 該当するノード・エッジが複数ある場合にリストとして返す
  • min/max/stddev/avg/sum: 数値処理
MATCH (p:Person)-[r]->(m:Movie)
WHERE  p.name = 'Tom Hanks'
RETURN m.title AS movie, type(r) AS relationshipType

exists

MATCH (p:Person)
WHERE exists ((p)-[:ACTED_IN]-())
SET p:Actor

split/coalesce

String型のプロパティ (e.g., "A|B") をリスト型 (e.g., ["A", "B"]) へ変換する事例が紹介されていた。

MATCH (m:Movie)
SET m.countries = split(coalesce(m.countries,""), "|"),
m.languages = split(coalesce(m.languages,""), "|"),
m.genres = split(coalesce(m.genres,""), "|")

collect

MATCH (p:Person)-[:ACTED_IN]->(m:Movie)
WHERE p.name ='Tom Cruise'
RETURN collect(m.title) AS tomCruiseMovies, size(collect(m)) AS tomCruiseMovieCount

shortestPath

MATCH p = shortestPath((p1:Person)-[*]-(p2:Person))
WHERE p1.name = "Eminem"
AND p2.name = "Charlton Heston"
RETURN  p
MATCH p = shortestPath((p1:Person)-[:ACTED_IN*]-(p2:Person))
WHERE p1.name = "Eminem"
AND p2.name = "Charlton Heston"
RETURN  p
MATCH (p:Person {name: 'Eminem'})-[:ACTED_IN*1..4]-(others:Person)
RETURN  others.name

APOC系

APOC (Awesome Procedures on Cypher) ライブラリの関数は厳密にはCypherの関数ではないので CALL で呼び出す必要がある

  • apoc.meta.nodeTypeProperties: 各ノードが有するプロパティの情報 (e.g., データ型) を返す
  • apoc.meta.relTypeProperties: 各エッジが有するプロパティの情報 (e.g., データ型) を返す
CALL apoc.meta.nodeTypeProperties()

Parameters

:paramでパラメータを設定、Cypherクエリの中で使うことが可能

// パラメータの設定
:param number: 10
:param actorName: 'Tom Hanks'

// 設定されたパラメータを確認する
:params: 
// :params {} // 全パラメータを削除

MATCH (p:Person)-[:ACTED_IN]->(m:Movie)
WHERE p.name = $actorName
RETURN p, m

JavascriptとかPythonとかで以下のように使える

def get_actors(tx, movieTitle): # (1)
  result = tx.run("""
    MATCH (p:Person)-[:ACTED_IN]->(m:Movie)
    WHERE m.title = $title
    RETURN p
  """, title=movieTitle)

  # Access the `p` value from each record
  return [ record["p"] for record in result ]

with driver.session() as session:
    result = session.read_transaction(get_actors, movieTitle="Toy Story")

Graph modeling

Naming conventions

ラベル・リレーション・プロパティには以下の命名規則を使うことがベストプラクティスとされている。

  • ラベル (Node type): CamelCase (e.g., Company, HighSchool)
  • リレーション (Edge type): 大文字 & アンダースコア (e.g., FOLLOWS, MARRIED_TO)
  • プロパティ: camelCase (e.g., deptId, firstName)

Importing CSV

LOAD CSV

  • LOAD CSV
  • WITH HEADERS: ヘッダー行が存在する場合に指定する
  • FIELDTERMINATOR: csvでない場合に使う。区切り文字を指定する。
CALL {
LOAD CSV WITH HEADERS
FROM 'https://data.neo4j.com/importing/2-movieData.csv'
AS row
WITH row WHERE row.Entity = "Join" AND row.Work = "Acting"
MATCH (p:Person {tmdbId: toInteger(row.tmdbId)})
MATCH (m:Movie {movieId: toInteger(row.movieId)})
MERGE (p)-[r:ACTED_IN]->(m)
ON CREATE
SET r.role = row.role
SET p:Actor
}

Neo4j Data Importer

GUIで動かせるツール。ローカル環境にあるCSV形式のデータ (1行目にヘッダが必須) をNeo4jへインポートできる。よくできているがエンジニアが使うかは微妙と感じる。リスト型のプロパティはインポートできないなどプログラムでのデータインポートと比較して制限がある。

*1:Undirected edgeの作り方がよくわからない。同じクエスチョンがGithubでissueとして投稿されている https://github.com/neo4j/neo4j/issues/13009

Pyroでベイズモデリング①:基本の確認

Pyro

PyroはUberが作った確率的プログラミングのためのPythonライブラリ。PyTorchをベースにしており、PyTorchで実装したニューラルネットを組み込んだ深層生成モデルの開発ができる。最近、シングルセルデータやマルチオミクスデータを対象とした深層生成モデルが多数報告されており、自分でも実装&既存のモデルをカスタマイズできるようになりたいと考えていたので勉強してみた。Pyroの公式ドキュメントは充実しており、読んでいるだけで面白く勉強になるのだが、試している中で疑問に思ったことが多々あったので、メモを残していこうと思う。

pyro.ai

Stanのような歴史の長い確率的プログラミング言語と比較すると、日本語で書かれている資料は少ない。多くの日本語の情報は「ベイズ線形回帰をPyroで実装してみた」みたいな内容が多いので、このシリーズでは、より複雑なモデルを実装することにトライしていこうと思う。

本記事では以下のページが参考になった。

実装の流れ

Step 1: ModelとGuideの定義

Pyroでは、はじめにModelとGuideを関数として実装する。Modelにはモデルと各パラメータの事前分布を、Guideには変分推論 (後述) で用いる近似事後分布を定義する。

Pyroのチュートリアルで紹介されている線形回帰の事例をもとに実装の流れをさらっていく。この事例では国土の地形の凹凸の多さ(TRI: Terrain Ruggedness Index)とGDPの関係を線形回帰で調べており、アフリカ以外の地域ではTRIが高い (地形の凹凸が高い) ほどGDPが低いのに対して、アフリカではそれと逆の傾向があることを示している。ここでは以下のモデルに対してベイズ推論を行う。変数 Africaはアフリカ地域であれば1が入るバイナリ変数である。

  •  log(GDP) \sim \rm{Normal}(\mu, \sigma)
    •  \mu = a + b_{A} * Africa + b_{R} * TRI + b_{AR} * Africa * TRI
      •  p(a) \sim \rm{Normal}(0, 10)
      •  p(b_{A}) \sim \rm{Normal}(0, 1)
    •  p(\sigma) \sim \rm{Uniform}(0, 10)

このケースであれば各パラメータの事後分布は解析的に求めることができるが、あえて変分推論で近似事後分布を求める。ここでは変分ベイズでよく用いられる平均場近似 (後述) という近似法を使う。事後分布 p(b_{A}, b_{R}, b_{AR} | X)を近似する、3つの独立した分布 ( q(b_{A}), q(b_{R}), q(b_{AR})) を探索する。以下のコードでは分布 q(b_{A}), q(b_{R}), q(b_{AR})にはそれぞれ正規分布を仮定している。

  •  p(b_{A}, b_{R}, b_{AR} | X) \approx q(b_{A}) q(b_{R}) q(b_{AR})

このときModelとGuide (近似事後分布) はPyroで以下のように実装できる。Model/Guideを構成する要素 (Primitives) については後述する。

import torch
from torch.distributions import constraints
import pyro
import pyro.distributions as dist

def model(is_cont_africa, ruggedness, log_gdp):
    a = pyro.sample("a", dist.Normal(0., 10.))
    b_a = pyro.sample("bA", dist.Normal(0., 1.))
    b_r = pyro.sample("bR", dist.Normal(0., 1.))
    b_ar = pyro.sample("bAR", dist.Normal(0., 1.))
    sigma = pyro.sample("sigma", dist.Uniform(0., 10.))
    mean = a + b_a * is_cont_africa + b_r * ruggedness + b_ar * is_cont_africa * ruggedness
    with pyro.plate("data", len(ruggedness)):
        pyro.sample("obs", dist.Normal(mean, sigma), obs=log_gdp)

def guide(is_cont_africa, ruggedness, log_gdp):
    a_loc = pyro.param('a_loc', torch.tensor(0.))
    a_scale = pyro.param('a_scale', torch.tensor(1.), constraint=constraints.positive)
    sigma_loc = pyro.param('sigma_loc', torch.tensor(1.), constraint=constraints.positive)
    weights_loc = pyro.param('weights_loc', torch.randn(3))
    weights_scale = pyro.param('weights_scale', torch.ones(3), constraint=constraints.positive)
    a = pyro.sample("a", dist.Normal(a_loc, a_scale))
    b_a = pyro.sample("bA", dist.Normal(weights_loc[0], weights_scale[0]))
    b_r = pyro.sample("bR", dist.Normal(weights_loc[1], weights_scale[1]))
    b_ar = pyro.sample("bAR", dist.Normal(weights_loc[2], weights_scale[2]))
    sigma = pyro.sample("sigma", dist.Normal(sigma_loc, torch.tensor(0.05)))

変分推論 (Variational inference)

変分推論については多くの書籍・ウェブサイトで説明されているので、言及する必要はないかもしれないが、自分の中で腑に落ちた説明を以下に残しておく。

N個のパラメータ  Theta = (\theta_{1},...,\theta_{N}) を含むモデルを考える。複雑なモデルである場合、パラメータの事後分布 p(\Theta|X)は解析的に求めることができないことが多い。*1 このとき変分推論ではパラメータ ( \Theta = (\theta_{1},...,\theta_{N})) の事後分布 p(\Theta|X)を近似する近似事後分布 q(\Theta)を探索する。以下のKLダイバージェンスが最小となる分布 q(\Theta)を探す。

 q(\Theta) = \rm{argmin}_{q(\Theta)} \rm{KL} \lbrack q(\Theta) || p(\Theta|X) \rbrack \tag{1}

ここでは平均場近似による近似がよく使われる。平均場近似では各パラメータ \theta_{i}の分布が互いに独立していて、以下が成り立つ条件を仮定する。

 q(\Theta) = \Pi_{i} q(\theta_{i}) \tag{2}

式(1)の最適解を探索するには以下の変分下限(ELBO)  L \lbrack q(\Theta) \rbrack を最大化する q(\Theta)を探索すれば良い。

 L \lbrack q(\Theta) \rbrack = \int q(\Theta) \ln \frac{p(X, \Theta)}{q(\Theta)} d\Theta \tag{3}

この根拠は以下のように説明できる。 p(X)は以下のように変形できる。左辺 \ln p(X)はデータ Xに対して一定の値なので、\rm{L} \lbrack q(\Theta) \rbrackを最大化することは、KLダイバージェンスを最小化することと等しいと考えられる。

{
\begin{align}
\ln p(X) &= \int q(\Theta) \ln p(X) d\Theta \tag{4} \\\
&= \int q(\Theta) \ln \frac{p(X, \Theta)}{p(\Theta|X)} d\Theta \\\
&= \int q(\Theta) \ln \frac{q(\Theta)p(X, \Theta)}{q(\Theta)p(\Theta|X)} d\Theta \\\
&= \int q(\Theta) \ln \frac{p(X, \Theta)}{q(\Theta)} d\Theta - \int q(\Theta) \ln \frac{p(\Theta|X)}{q(\Theta)} d\Theta \\\
&= \rm{L} \lbrack q(\Theta) \rbrack + \rm{KL} \lbrack q(\Theta) || p(\Theta|X) \rbrack
\end{align}
}

Primitives

Pyroでのモデルの構築に使われる基本単位 (Primitivesと呼ばれる) についてまとめる。

pyro.param

  • パラメータを定義する際に用いる関数
  • 引数
    • name: パラメータの名前 (ParamStoreDictでの管理に使われる)
    • init_tensor: 初期値
    • constraint: 制約条件 (e.g. 正値)
  • パラメータはグローバルで管理されParamStoreDictと呼ばれる領域に保存される。*2
    • pyro.get_param_store()ParamStoreDictの内容をdict型で取得できる
    • pyro.clear_param_store()ParamStoreDictの内容を消去できる。学習後のパラメータを消去したい際に用いる。
  • ParamStoreDictのkeyにはpyro.paramの第一引数"name"の情報が使われるので、パラメータを宣言する際は、他のパラメータと同じ"name"を使わないようにしなければならない。
a_loc = pyro.param('a_loc', torch.tensor(0.))
d = pyro.get_param_store()
d["a_loc"]
tensor(0., requires_grad=True)

パラメータは以下の方法でも取得できる。

pyro.param("a_loc").item()

pyro.sample

  • 確率変数を定義する際に用いる関数
  • 引数
    • name: 確率変数の名前。これをもとにModelの分布とGuideの近似分布を対応させる
    • fn: 分布 (pyro.distributionsのクラスが用いられる)
    • obs: 実測値 ELBOの計算に使われる (default: None)
  • 返値
    • obsが空(None)の場合: 分布fnからサンプリングした値を返す
    • obsが空でない場合: obsに指定した値が返される
a = pyro.sample("a", dist.Normal(0., 1.))
a
tensor(0.5161)

pyro.plate

  • 繰り返し(グラフィカルモデルにおけるプレート)を定義する関数
    • e.g. N個の遺伝子の発現量を1つずつ同じ分布でモデル化する場合
    • pyro.plate内部の確率変数は(特定の条件下で)独立とみなされる
      • 例えば以下のコードでは特定のmean/sigmaにおいてobs内の各変数は独立といえる
      • 並列化により計算を高速化できる
  • 引数
    • name: プレートの名前
    • size: 繰り返しのサイズ
    • subsample_size: ミニバッチのサイズ (default: size)
    • dim: plateをどの次元に対応させるか (default: -1)
      • dim=-1の場合はbatch_shape (後述) のうち一番後ろ (後ろから1番目) をPlateに対応させる
mean = pyro.sample("mean", dist.Normal(0., 1.))
sigma = pyro.sample("sigma", dist.Uniform(0., 10.))

with pyro.plate("data", 5) as idx:
    y = pyro.sample("y", dist.Normal(mean, sigma))
y.shape
torch.Size([5])

batch_shapeとevent_shape

pyro.distributionモジュールの分布クラス (e.g. dist.Normal) は、パラメータに複数の値を与えたり、メソッドexpandを使うことで、複数の分布を同時に定義することができる。このとき、メソッドsampleや、上述のpyro.sampleにより、複数の分布からサンプリングされた値をtorch.tensor型のオブジェクトとして取得することができる。例えば、正規分布 (dist.Normal) のパラメータ ( \mu, \sigma) に4つの値を与えてると、4つの正規分布が同時に定義され、各分布からサンプリングされた値をサイズ4のテンソルで取得することができる。

means = torch.tensor([0.0, 1.0, 2.0, 3.0])
sds = torch.tensor([1.0, 1.0, 1.0, 1.0])
d = dist.Normal(means, sds)
x = pyro.sample("x", d)
x
tensor([0.2222, 1.5773, 0.2363, 3.2264])

これに関連する重要な概念として、分布クラスはbatch_shapeevent_shapeという属性を持つ。batch_shapeは分布の数を表し、event_shapeは各分布から出力される値の次元数を表す。分布クラスからサンプリングされた値の次元は「batch_shape + event_shape」の形になる (event_shapeが末尾の次元に対応する)。

以下の2つのケースは似ているが、異なるbatch_shapeevent_shapeとなる。①4つの正規分布を作る場合は、batch_shapeが[4]に対して、event_shapeは (=[1])となる。②の4次元正規分布を作る場合は、①と同様に4次元の値を返す分布が作られるが、1つの分布として扱われるため、batch_shapeが、event_shapeが[4]となる。

① 4つの正規分布を作る場合

means = torch.tensor([0.0, 1.0, 2.0, 3.0])
sds = torch.tensor([1.0, 1.0, 1.0, 1.0])
d = dist.Normal(means, sds)
d.batch_shape, d.event_shape, d.sample().shape
(torch.Size([4]), torch.Size([]), torch.Size([4]))

② 1つの4次元正規分布を作る場合

means = torch.tensor([0.0, 1.0, 2.0, 3.0])
sds = torch.eye(4)
d = dist.MultivariateNormal(means, sds)
d.batch_shape, d.event_shape, d.sample().shape
(torch.Size([]), torch.Size([4]), torch.Size([4]))

to_event

分布クラスのメソッドto_eventを使うと、複数の独立した確率分布を、一つの多次元確率分布として扱うことができるようになる。

① 3x4の独立したベルヌーイ分布

d = dist.Bernoulli(0.5 * torch.ones(3,4))
d.batch_shape, d.event_shape, d.sample().shape
(torch.Size([3, 4]), torch.Size([]), torch.Size([3, 4]))

② 3つの独立した4次元ベルヌーイ分布

d = dist.Bernoulli(0.5 * torch.ones(3,4)).to_event(1) # 末尾1次元を非独立として扱う
d.batch_shape, d.event_shape, d.sample().shape
(torch.Size([3]), torch.Size([4]), torch.Size([3, 4]))

pyro.plateの引数dimについて

上述の通りpyro.plateの引数dimはplateがどのbatch_shapeに対応するかを指定するのに使われる。例えば、以下のコードでは、event_shapeが[3, 4]のベルヌーイ分布が5つ (batch_shape=5) 作られる。サンプリングされる値の次元は「batch_shape + event_shape」の形になるので、この場合は[5, 3, 4]次元のテンソルが出力される。

with pyro.plate("a", 5, dim=-1):
    d = pyro.sample("d", dist.Bernoulli(0.5 * torch.ones(3,4)).to_event(2))
d.shape
torch.Size([5, 3, 4])

Step 2: 学習

ModelとGuideを実装したあとはパラメータの推論を行う。Pyroの変分推論では、Adamなどの確率的な最適化アルゴリズムを使って、事後分布にできるだけ近い近似分布を探索する (SVI: Stochasticなvariational inference)。はじめに、pyro.optimモジュールのクラスを使い最適化のアルゴリズム・条件 (e.g. learning rate) を設定し、SVIクラスでModel/Guide・Optimizer・損失関数 (この場合はELBO) を指定する。つづいてsvi.stepでパラメータの更新を行う。svi.stepはModel/Guideの引数となるデータを受け取り、パラメータを更新し、返り値として損失 (ELBO) を出力する。以下の事例では、for文で5000回パラメータの更新を行っている。

from pyro.infer import SVI, Trace_ELBO
from pyro.optim import Adam

pyro.clear_param_store() # パラメータのリセット

optimizer = Adam({"lr": .05})
svi = SVI(model, guide, optimizer, loss=Trace_ELBO())

num_iters = 5000
for i in range(num_iters):
    elbo = svi.step(is_cont_africa, ruggedness, log_gdp)
    if i % 500 == 0:
        print("Elbo loss: {}".format(elbo))

Step 3: 事後分布の取得

pyro.infer.Predictiveを使うと以下のように事後分布から各パラメータをサンプリングできる。

from pyro.infer import Predictive
params = ["a", "bA", "bR", "bAR"]
pred_model = Predictive(model, guide, num_samples=2000, return_sites=params)
pred_sample = pred_model(is_cont_africa, ruggedness, log_gdp)
pred_sample["a"]
tensor([[9.1086],
        [9.1260],
        [9.1357],
        ...,
        [9.2670],
        [9.0693],
        [9.0717]])

*1:https://ntacoffee.com/variational-inference/

*2:どうしてこのようなデザインになったのだろうか

シングルセル解析⑧: 細胞の運命を確率的に推定する (Plantir, CellRank)

問題設計

CellRankは、Fabian TheisとDana Pe'erのグループから最近報告された、シングルセルデータから細胞が将来的に辿る状態遷移の可能性を確率的に求める手法。「細胞Aは細胞Bよりも細胞Cへと辿り着く可能性が高い」といった推定ができる。

www.nature.com

Plantir

CellRankのアイデアはDana Pe'erのグループから2019年に発表されたPlantirから発展したものだと思われる。Plantirは、細胞の状態遷移にマルコフ性がある (現在の状態のみが次の状態に影響し、過去の状態は影響しない) と仮定し、上記の課題に挑んでいる。手法の概要は以下の通り*1

  1. Diffusion mapでscRNA-seqデータを次元圧縮
  2. Top nのDiffusion components (DCs) を使いkNN graph  G_{E}を求める
  3. 各細胞のPseudo-time  \tau (解析者が選んだ起点からのkNN graphにおける最小パス数) を求める
  4. Pseudo-timeをもとに G_{E}を有向グラフ G_{D}に変換する
    •  \tau_{i} - \tau_{j} \leq \sigma_{i}の場合は細胞 iから細胞 jへのエッジと考える*2
  5. 解析者が選んだ終点からのEdgeを G_{D}から除外 *3
  6.  G_{D}から細胞間の状態遷移確率行列 Pを求める
  7. 状態遷移確率行列 Pを使ったランダムウォーク

結果として下図dのような結果が得られる。例えば、図dのTrajectoryの(2)に位置する細胞は細胞状態A/B/Cのいずれにも同程度分化できるポテンシャルがあるが、(6)の細胞はB/Cにのみ分化する、といったような解釈ができる。

引用元:Setty et al. “Palantir characterizes cell fate continuities in human hematopoiesis” bioRxiv 385328; doi: https://doi.org/10.1101/385328. CC-BY 4.0

www.nature.com

CellRank

CellRankはマルコフ連鎖での細胞状態遷移のモデリングRNA velocity*4を利用する。RNA velocityを使うことによって、(i) 解析者がTrajectoryの始点・終点を決める必要がなく、(ii) Pseudo-timeに依存したPlantirのように「始点から終点へ」という単一方向の状態遷移を仮定する必要がなくなる。したがって、より複雑な細胞状態遷移の過程を調べるのに適した手法だと言える。RNA velocityに関しては以前の記事を参考のこと。手法の概要は以下の通り。7-8のあたりの展開は正直理解できていない。

  1. kNN graphを作る
    • Plantirと違いDiffusion map以外の次元圧縮法を適用可能 (Default: PCA)
  2. Cell-cell adjacency matrix  A \in R^{N \times N}
    • 細胞 i jがいずれかのkNNに含まれていない場合:  A_{ij} = 0
    • それ以外:  A_{ij} = \rm{dist}(x_{i}, x_{j})
      • dist: デフォルトはユークリッド距離
      •  x_{i}: first top n PCs of cell  i (Default: top 30 PCs)
  3. Connectivity kernel  P^{s} \in R^{N \times N}
    •  P^{s} = row-normalized  A
  4. Velocity kernel  P^{v} \in R^{N \times N}
    • 下図b参照
    •  P^{v}_{ij} = \frac{ \rm{exp} (\sigma c_{ij}) }{ \sum_{k \in K_{i}} \rm{exp} (\sigma c_{ik}) }
      •  c_{ij} = \rm{corr} (v_{i}, s_{ij})
        •  v_{i}: 細胞 iのVelocity
        •  s_{ij} = x_{i} - x_{j}: 細胞 i jの遺伝子発現の差
      •  K_{i}: 細胞 iのkNN
      •  \sigma = 1 / \rm{median}_{k \in K_{i}}(c_{ik})
        • VelocityがkNNのどの方向にも向いていない場合に P^{v}_{ij}を小さくする用途のパラメータ
  5. Merged kernel  P \in R^{N \times N}
    •  P = (1 - \lambda) P^{v} + \lambda P^{s}
      • 論文では \lambda=0.2が安定した結果を出すと記載されている
      • カーネルの設定はユーザーに委ねられている*5。ベストプラクティスが特に無いので \lambda=0.2を使うしかなさそう。
  6. Macrostate identification
    • 単一細胞レベルのデータであるカーネル Pはノイズが大きいので、細胞集団(Macrostateと呼ぶ)レベルの情報へと変換する
    • GPCCA (Generalized Perron Cluster Cluster Analysis)
      1. カーネル PSchur分解
        •  P = QRQ^{T}
        •  Q \in R^{N \times N}: ユニタリ行列 ( QQ^{*} = I) 各列がMacrostateに対応する基底ベクトル
        •  R \in R^{N \times N}:  P固有値 \lambda_{1}...\lambda_{l}が対角成分にある上三角行列
      2. 以下のいずれかの基準で上位のMacrostateを選ぶ
        • Eigengap heuristics (固有値閾値を設ける方法) *6
        • minChi algorithm
        • Crispness
      3. 選ばれたMacrostateのみを含む基底ベクトル行列を \hat{Q} \in R^{N \times n_{s}}と呼ぶ
  7. Macrostate membership  \chi \in R^{N \times n_{s}}
    •  \chi = \hat{Q}A:  \chi_{is}は細胞 iのState  sとの類似度を表す
      •  A: Rotation matrix *7
    •  Aは以下の f_{n_{s}}(A)を最小化する最適化問題を解き求める
      •  f_{n_{s}}(A) = n_{s} - \rm{trace}( \hat{D}^{-1} \chi^{T} D \chi)
      •  \hat{Q}^{T} D Q = I
      • 制約条件
        •  \chi_{is} \leq 0
        •  \sum_{s} \chi_{is} = 1
      • これにより \chiの各列 (各Macrostateに該当する細胞のパターン) は直交に近づく
  8. Macrostate-level transition matrix  P^{c} \in R^{n_{s} \times n_{s}}
    •  P^{c} = ( \chi^{T} D \chi )^{-1} (\chi^{T} D P \chi)
  9. Terminal stateを推定する
    • 分化の終末にいるMacrostateは他の細胞への状態遷移が少ないと考えられる。したがってCellRankでは遷移行列 P^{c}の情報をもとに終末にあるMacrostateを以下のように定める。
    • State  i is terminal if  P^{c}_{ii} > \epsilon
      •  \epsilon: Threshold (Default: 0.96)
  10. Initial stateを推定する
    1. Stationary distribution  \pi \in R^{N}を求める
      •  \pi: 定常状態にある細胞状態組成
      •  \pi^{T}P = \pi^{T}
      •  \sum_{i} \pi^{T}_{i} = 1
    2.  \piをMacrostate-levelに変換
      •  \pi_{p} = \chi^{T} \pi
    3. Initial state
      • 定常状態を「一定以上の時間が経過した後の細胞状態組成」と仮定すると、Initial stateの細胞は定常状態では少なくなると考えられる。したがってCellRankでは以下の基準で選択する。Initial stateの数 nはユーザーが決める必要がある。
      • State  i is initial if  \pi_{i} is in top  n smallest values
  11. Terminal cells  R_{t}を定義
    • Terminal state  t ごとにMembership  \chi_{t}の値が高い細胞をTerminal cells  R_{t}とする
  12. 各細胞の運命を推定する
    • Transient cellとRecurrent cellという2つの状態を考える。Transient cellからはRecurrent cellへ遷移できるのに対し、反対の遷移は生じないとすると、Transition matrix  Pは以下のように表せる
      • このとき A=(I - Q)^{-1}S A_{ij}がTransient cell  iがRecurrent cellのうち jへ最初に到着する確率を表すことが証明されているらしい*8
    • CellRankではTerminal cell R_{t}をRecurrent cellと見なし、その他の細胞(Transient cell)への遷移をゼロに変換。上述の Aを求め、各細胞がどのTerminal stateへ辿り着くか (Absorption probabilityと呼んでいる)を推定する
 
P = \begin{bmatrix}
\hat{P}&0\\
S&Q
\end{bmatrix}

引用元:Lange et al. “ CellRank for directed single-cell fate mapping” bioRxiv 2020.10.19.345983; doi: https://doi.org/10.1101/2020.10.19.345983 CC-BY 4.0

*1:簡略化しているので詳細は論文を要参照

*2:σは

*3:ランダムウォークの際に終点で収束するように

*4:過去記事: https://auroratummy.hatenablog.com/entry/2021/12/07/011122

*5:https://cellrank.readthedocs.io/en/stable/kernels_and_estimators.html#Initialize-a-kernel

*6:https://cellrank.readthedocs.io/en/stable/kernels_and_estimators.html#Compute-a-matrix-decomposition

*7:前述のAdjacency matrixと違うので注意

*8:Tolverらの教科書が引用されていた。http://web.math.ku.dk/noter/filer/stoknoter.pdf

Context-specific metabolic modelingについて

Context-specific metabolic modeling

Context-specific metabolic modelingとは、トランスクリプトーム・プロテオーム等の発現情報を用いて、Recon3DやHuman1等の全代謝反応を含むGenome-scale metabolic model (GEMs)から、解析対象 (e.g. 特定の細胞、組織) で機能している・していない代謝反応のみを含むGEMsを抽出する手法。

2008年頃から色々な手法が提案されており、特に新しい研究分野ではないのだが、個人的なプロジェクトで似たような方法を試してみたいと考えている。ほとんどがMATLABで実装されていること (筆者はオープンソースで実装したい)、手法に応じて結果が大きく異なることが報告されていること *1 など、手法を理解していない状態で試すと危ない雰囲気なので、簡単にまとめてみようと思う。

Context-specific metabolic modelingの手法は以下のレビュー論文では3つのグループに分類されている。ここでは各グループの代表的な手法について調べる。

doi.org

GIMME-like family

上記論文では「発現量の低い遺伝子・タンパク質に対応する代謝反応の流束を最小化する」制約条件を利用した手法がGIMME-likeと分類されている。

GIMME

GIMMEは2008年にPalssonらによって発表された最初期の手法。(i) 遺伝子発現 (or タンパク発現) データ、(ii) 対象生物種のGEM、(iii) 発現の有無を決めるための閾値 (e.g. この閾値を上回る遺伝子は「発現あり」、下回る場合は「発現なし」とみなされる)、(iv) Required Metabolic Functionalities (RMF) (生物学的な知識に基づき存在していると考えられる代謝機能) を入力として受け取る。原著論文では「培地中の炭素源が生育につながる*2」ことをRMFとして設定している。

以下のステップで組織・細胞特異的なGEMsを構築する。

  1. 発現の無い (閾値を下回る) 遺伝子・タンパク質に対応する代謝反応をGEMから除外する
  2. RMFを実現できるように除外した代謝反応の一部を元に戻す
  3. RMFを最大化する代謝流束をFBAで推定する
  4. RMFがある値以上 (論文では3.で求めた最大値の0.9倍以上を基準としている) となる範囲で以下のInconsistency Score (IS) が最小となる代謝流束をFBAで計算する

Inconsistency scoreは下図のように設定されている。赤い矢印は対応する遺伝子の発現量が閾値 *3 を下回っている代謝反応を表す。最適な閾値を決める方法は本論文では特に提案されていない。閾値をいくつか試して結果への影響が小さかったことを確認している、とのこと。*4 上記のようにGIMME-like familyは遺伝子・タンパク質発現とある程度齟齬のない代謝流束で、Biomass functionなどの目的関数を最大化することを目指す手法である。

  •  \sum_{i} c_{i} v_{i}: Inconsistency Score
  •  v_{i}: 反応 iの流束
  •  c_{i} = max(0, x_{cutoff} - x_{i})
    •  x_{cutoff}: 閾値
    •  x_{i}: 反応 iに対応する遺伝子の発現量

引用元:Becker SA, Palsson BO. Context-specific metabolic networks are consistent with experiments. PLoS Comput Biol. 2008 May 16;4(5):e1000082. doi: 10.1371/journal.pcbi.1000082. PMID: 18483554; PMCID: PMC2366062.

doi.org

iMAT-like family

上記レビュー論文では「低発現遺伝子に対応する代謝反応をできるだけ削り、高発現遺伝子に対応する代謝反応をできるだけ残す」最適化を行う手法をiMAT-likeと呼んでいた。GIMME-like familyは遺伝子発現データとの齟齬の少ない制約条件下でBiomass functionなどのRMFを最大化する代謝状態を探索する方法であったが、iMAT-like familyの手法はRMFを考慮せず、シンプルに遺伝子発現データとの齟齬の少ない代謝状態を探索する。

iMAT

iMATは2008年にShlomi・Palsson・Ruppinらによって提案された手法。(i) 遺伝子発現 (or タンパク発現) データ、(ii) 対象生物種のGEM、(iii) 発現レベルをHigh/Normal/Lowに分けるための閾値、(iv)パラメータ \epsilon (後述)を入力として受け取る。

以下の制約条件下で目的関数が最大となる代謝流束を見つける。この目的関数を最大化することで、高発現遺伝子に対応する代謝の流束ができるだけ非ゼロの値になり、低発現遺伝子に対応する代謝の流束がゼロになるような、遺伝子発現データとの齟齬の少ない代謝状態を探索できる。下記論文のFig. 1がイメージをつかみやすい*5。この問題はパラメータ yが整数値なので混合整数線形計画法 (MILP)で最適化を行う。MILPは線形計画法(LP)よりも最適解の探索が困難であるため*6、おそらくGIMMEなどのLPベースの手法よりも計算時間を要すると思われる。

  • 目的関数:  \sum_{i \in R_{H}} (y_{i}^{+} + y_{i}^{-}) + \sum_{i \in R_{L}} \hat{y}_{i}^{+}
    •  R_{H}: 高発現遺伝子に対応する代謝反応
    •  R_{L}: 低発現遺伝子に対応する代謝反応
    •  y_{i}^{+}: 代謝 iの順方向のFluxがある場合1、ない場合0のバイナリ変数
    •  y_{i}^{-}: 代謝 iの逆方向のFluxがある場合1、ない場合0のバイナリ変数
    •  \hat{y}_{i}^{+}: 代謝 iのFluxがある場合0、ない場合1のバイナリ変数
  • 制約条件
    1.  Sv = 0 *7
    2.  v_{min,i} \leq v_{i} \le v_{max,i} *8
    3.  v_{i} + y_{i}^{+} (v_{min,i} - \epsilon) \geq v_{min, i}, i \in R_{H}
    4.  v_{i} + y_{i}^{+} (v_{max,i} + \epsilon) \leq v_{max,i}, i \in R_{H}
      •  \epsilon: このパラメータを設けることで、代謝が存在する ( y_{i}^{+} or y_{i}^{+}= 1の場合)にFluxがある程度の値になるように制約をかける
    5.  v_{min,i}(1-\hat{y}_{i}^{+}) \leq v_{i} \leq v_{max,i}(1-\hat{y}_{i}^{+}), i \in R_{L}
    6.  y_{i}^{+}, y_{i}^{-}, \hat{y}_{i}^{+} \in [0, 1]

www.nature.com

INIT (tINIT)

INITはJens Nielsenらのグループが2012年に発表した手法。その後2014年に改良版のtINITが発表された。INITは以下のように定式化されている。原著論文を読んでいてもわからない点があったのでレビューを参考にした。iMATと同様に一部のパラメータが整数値なのでMILPで最適化を行う。

  • 目的関数:  \sum_{i \in R} w_{i} y_{i}
    •  R: Reactions
    •  M: Metabolites
    •  w_{i}: 代謝反応 iのWeight
      • 代謝反応 iに対応する遺伝子・タンパクの発現レベルに応じて値を変える
      • 2012年のオリジナル論文では w_{i} = 5log(\frac{Average_{i}}{Signal_{i,j}})としている *9
  • 制約条件
    •  Sv = b: INITは定常状態を仮定しない
      •  b_{j} \geq \sigma, j \in KeyMet: 代謝 jが対象細胞・組織に存在することが知られている場合は \sigma *10 以上物質産生が起こっていると仮定 *11
      •  b_{j} = 0, j \notin KeyMet
    •  y_{i} \in [0,1]: 代謝反応  i が存在する場合は1、除外された場合は0
    •  y_{i} v_{min,i} \leq v_{i} \leq y_{i} v_{max,i}

引用元:Agren R, Bordel S et al. Reconstruction of genome-scale active metabolic networks for 69 human cell types and 16 cancer types using INIT. PLoS Comput Biol. 2012;8(5):e1002518. doi: 10.1371/journal.pcbi.1002518. Epub 2012 May 17. PMID: 22615553; PMCID: PMC3355067.

journals.plos.org

https://www.embopress.org/doi/full/10.1002/msb.145122www.embopress.org

MBA-like family

遺伝子発現データを元にGEMから発現の無い枝を取り除く(Pruning)する手法。この手法はGIMME-like familyやiMAT-like familyのように特定の代謝状態を探索するのではなく、

MBA

MBAは2010年にJerby・Shlomi・Ruppinによって提案された手法。下記の目的関数を最大化するPartial model (PM)を以下のアルゴリズムで探索する*12。Globam model (GM) から少しずつ不必要な反応を削っていくアプローチである。2. でCore reactions ( C_{H}, C_{M}) に含まれない反応の集合  Pを抽出。3.以降のFor loopでは、 Pに含まれる反応 rを一つずつKOして、その結果Inactiveになる反応を除外する (枝葉を切る, pruning)。CheckModelConsistencyはモデル R_{P}から反応 rを除外したときにInactiveになる反応を返す関数。

2.において Pに含まれる反応 rの順番はランダムに決める。反応 rの並びによって最終的なGEMの形が変わるので、論文ではpruningを1000回行い、できた1000個のGEMから共通性の高い領域を抽出している。

  • 目的関数:  |R_{P} \cap C_{M}| - \epsilon | R_{P} \setminus (C_{M} \cup C_{H})|
    • Core reactions:  C
      •  C_{H}: 高確率で存在する代謝反応 (e.g. 文献・実験由来の情報)
      •  C_{M}: 中確率で存在する代謝反応 (e.g. オミクス由来の情報)
    • Global model:  GM = (M_{G}, R_{G}, S_{G}, L_{G})
      • M: 代謝
      • R: 代謝反応
      • S: Stoichiometry matrix
      • L: Limitation of flux (各流束の最大・最小値)
    • Partial model:  PM = (M_{P}, R_{P}, S_{P}, L_{P})
      •  C_{H} \in R_{P}: PMは C_{H}の全反応を含む

アルゴリズム

f:id:auroratummy:20220118002936p:plain 引用元:Jerby, Livnat et al. “Computational reconstruction of tissue-specific metabolic models: application to human liver metabolism.” Molecular systems biology vol. 6 (2010): 401. doi:10.1038/msb.2010.56

www.ncbi.nlm.nih.gov

mCADRE

mCADREはMBAを改良した方法で2012年に発表された。MBAがFor loopのiterationの順番をランダムに決めていたことを問題視。近傍にCore reactionがあるかを考慮したスコアを設定し、Non-core reaction  Pに含まれる反応を並べることで、一意な結果が得られる(かつMBA論文のように1000回 iterationを実施する必要がないのでより短い時間で計算できる)手法を提案。

引用元:Wang, Y., Eddy, J.A. & Price, N.D. Reconstruction of genome-scale metabolic models for 126 human tissues using mCADRE. BMC Syst Biol 6, 153 (2012). https://doi.org/10.1186/1752-0509-6-153

bmcsystbiol.biomedcentral.com

Pythonで上記の手法を試すには

Pythonでは以下のライブラリに各手法が実装されている。Context-specific modelingの関数・クラスについての説明がほとんど存在しないツールもあるので、本当に使えるのか?だが、実装する際に参考になると思われる。

*1:https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1006867

*2:Biomass産生につながる

*3:この論文の場合は12.0

*4:Data not shownだが

*5:著作権の関係で引用できない

*6:https://www.msi.co.jp/nuopt/introduction/algorithm/mip.html

*7:一般的なFVAの制約条件:定常状態

*8:これも一般的なFVAの制約条件:事前情報に基づく各代謝の最大・最小流束

*9:Tissue jにおける発現値を全Tissueの平均値と比較

*10:small positive value

*11:これによってこの物質の産生に関わる代謝反応が除外されることを防ぐ

*12:なぜかこの論文では差集合をバックスラッシュでなくスラッシュで表しているので注意

マイクロバイオーム解析④:COBRApyで細菌の代謝をモデリングする(2)

2細菌が同じ環境にいるときの代謝を推定する

前回の記事の続き。2種類の細菌が同じ環境にいるときに、Cooperative/Competitiveな関係になるのか、どのように代謝物を交換するのか、といった細菌間相互作用をCOBRApyを使ったFBA (Flux Balance Analysis)で推定してみる。

ここでは以下の論文のFig. 4dで示されているように、Bacteroides caccae ATCC 43185とLactobacillus rhamnosus GGをDMEM 6429という培地で共培養した場合の細菌間相互作用を推定してみる。

www.nature.com

本記事の検証に使ったコードはgithubにアップした。

github.com

2生物種のGEMsを作る方法

この論文ではBacteroides caccae ATCC 43185とLactobacillus rhamnosus GGのGEMs (Genome-scale metabolic models) を連結させた代謝モデルを作り、両細菌の生育率の合計が最大になる場合の代謝流束をFBAで求めている。メソッドを読むと筆者らが2015年の論文に用いた細菌間に共通する区画*1を設ける手法を利用していた。この手法のオリジナルは以下の2010年の論文のようだ。このオリジナル論文のFig. 2が仕様を理解する上でわかりやすかった。

doi.org

引用元:Kitgord and Segre. PLoS Comput Biol. 2010 Nov 18;6(11):e1001002. doi: 10.1371/journal.pcbi.1001002.

Fig. ABは各細菌のGEMsである。区画CYT1CYT2は各細菌の細胞内領域、EX1EX2は各細菌の細胞外領域、XYZ代謝物質を表す。この手法ではFig. Cのように2つのGEMsに追加して、細菌間で共通の区画ENVを設定する。

これはStoichiometric matrixで表すとFig. EFGのように表せる。各細菌のStoichiometric matrix (Fig. EF) に若干手を加えるだけで2種共培養のStoichiometric matrix (Fig. G) を作ることができる。

MATLAB版のCOBRAには2細菌のモデルを作るための機能が用意されているが*2、COBRApyでは見当たらない。すでにMMINTEというライブラリが同様の機能をPythonで実装しているようなので、それを使ってみる。

doi.org

MMINTEのgithubページにおいてあったスクリプトを参考に以下のコードを走らせるとエラーが出た。

import mminte
species = ["Bacteroides_caccae_ATCC_43185", "Lactobacillus_rhamnosus_GG_ATCC_53103"]
model_filenames = ["AGORA/CurrentVersion/AGORA_1_03/AGORA_1_03_sbml/{s}.xml".format(s=s) for s in species]
pairs = mminte.get_all_pairs(model_filenames)
pair_model_filenames = mminte.create_interaction_models(pairs, ".")
pair_model_filenames

RuntimeError: The heuristic for discovering an external compartment relies on names and boundary reactions. Yet, there are neither compartments with recognized names nor boundary reactions in the model.

エラーメッセージを見ても原因がよくわからなかったので、MMINTEのコードをチェックしてみた。上記の関数mminte.create_interaction_modelsの中では関数create_community_modelが主要な働きをしている様子。

import six
community_compartment = "u" # lumen

# MetaboliteのSuffixをチェックする
#   e.g. ****[e] はExtracellular metaboliteを表すbigg形式のSuffix
modelseed_suffix_re = re.compile(r'_([ce])$')
bigg_suffix_re = re.compile(r'\[([cepu])\]$')
def _id_type(object_id):
    if re.search(bigg_suffix_re, object_id) is not None:
        return 'bigg'
    if re.search(modelseed_suffix_re, object_id) is not None:
        return 'modelseed'
    return 'unknown'

# Metaboliteの場所 (Compartment) を変更する関数
def _change_compartment(id_str, compartment, id_type):
    if id_type == 'bigg':
        return re.sub(bigg_suffix_re, '[{0}]'.format(compartment), id_str)
    elif id_type == 'modelseed':
        return re.sub(modelseed_suffix_re, '_{0}'.format(compartment), id_str)
    return '{0}_{1}'.format(id_str, compartment)


def create_community_model(model_files, model_names):
    # 空のGEMsを作成する
    community = cobra.Model('Community')

    # Exchange reaction IDのリスト (後に重複しないようにリストを作成する)
    exchange_reaction_ids = set()

    # 各生物のGEMsを追加する
    for model_file, model_name in zip(model_files, model_names):
        # GEMの読み込み
        model = cobra.io.read_sbml_model(model_file)

        # 
        linear_coefficients = cobra.util.solver.linear_reaction_coefficients(model)
        objective_id = '{0}_{1}'.format(model_name, six.next(six.iterkeys(linear_coefficients)).id)
        objective_value = six.next(six.itervalues(linear_coefficients))

        # 各MetaboliteのCompartmentの情報が存在するかMetabolite IDのSuffixをもとに確認
        for metabolite in model.metabolites: metabolite.notes['type'] = _id_type(metabolite.id)
        unknown = model.metabolites.query(lambda x: 'unknown' in x['type'], 'notes')
        if len(unknown) > 0:
            raise Exception('Unknown compartment suffixes found in metabolites for {0}'.format(model_file))

        # 各細菌のGEMsのExchange reactionのリストを取得する
        exchange_reactions = model.exchanges

        # 細菌間で共通の区画`ENV`を設定する
        for reaction in exchange_reactions:
            if reaction.id not in exchange_reaction_ids:
                exchange_reaction_ids.add(reaction.id)
                metabolite = six.next(six.iterkeys(reaction.metabolites)).copy()
                # MetaboliteのCompartmentを"e"からlumen ("u")へ変更
                metabolite.compartment = community_compartment
                metabolite.id = _change_compartment(metabolite.id, community_compartment, metabolite.notes['type'])
                community.add_metabolites(metabolite)
                rxn = community.add_boundary(metabolite)  # This is slow on cobra06
                rxn.id = reaction.id  # Keep same ID as species model

        for reaction in model.reactions:
            if reaction in exchange_reactions:
                # 区画`ENV`と`EX`の間の物質交換を定義する
                e_metabolite = six.next(six.iterkeys(reaction.metabolites))
                u_metabolite = community.metabolites.get_by_id(
                    _change_compartment(e_metabolite.id, community_compartment, e_metabolite.notes['type']))
                reaction.add_metabolites({u_metabolite: 1.})
                # Reaction ID の名前を変える (各細菌の反応を区別するため)
                reaction.id = '{0}_TR_{1}'.format(model_name, reaction.id[3:]) # Strip "EX_" prefix
                reaction.bounds = (-1000., 1000.)
            else:
                # Reaction ID の名前を変える (各細菌の反応を区別するため)
                reaction.id = '{0}_{1}'.format(model_name, reaction.id)

        # Metabolite ID の名前を変える (各細菌の物質を区別するため)
        for metabolite in model.metabolites:
            if metabolite.compartment != community_compartment:
                metabolite.id = '{0}_{1}'.format(model_name, metabolite.id)

        # あたらしいGEMに追加する
        community += model

        # 各細菌のBiomass functionを目的関数に追加
        objective_reaction = community.reactions.get_by_id(objective_id)
        cobra.util.solver.set_objective(community, {objective_reaction: objective_value}, additive=True)
    
    # GEMに名前をつける
    community.id = 'x'.join(model_names)

    return community

cobra.Model.add_metabolitesで上記のエラーメッセージが出ていた。調べると以下のメソッドが関わっていそうなことがわかった。cobra.Model.compartments == "u"が細胞外を表す区画だと認識されていないようだ。おそらくMMINTEが開発された頃からCOBRAの仕様が変わったものと思われる。

import logging
LOGGER = logging.getLogger(__name__)

# cobra.Model.compartmentsから細胞外のCompartmentを見つける
def find_external_compartment(model):
    if model.boundary:
        counts = pd.Series(tuple(r.compartments)[0] for r in model.boundary)
        most = counts.value_counts()
        most = most.index[most == most.max()].to_series()
    else:
        most = None
    like_external = cobra.medium.annotations.compartment_shortlist["e"] + ["e"]
    matches = pd.Series(
        [co in like_external for co in model.compartments], index=model.compartments
    )

    if matches.sum() == 1:
        compartment = matches.index[matches][0]
        LOGGER.info(
            "Compartment `%s` sounds like an external compartment. "
            "Using this one without counting boundary reactions" % compartment
        )
        return compartment
    elif most is not None and matches.sum() > 1 and matches[most].sum() == 1:
        compartment = most[matches[most]][0]
        LOGGER.warning(
            "There are several compartments that look like an "
            "external compartment but `%s` has the most boundary "
            "reactions, so using that as the external "
            "compartment." % compartment
        )
        return compartment
    elif matches.sum() > 1:
        raise RuntimeError(
            "There are several compartments (%s) that look "
            "like external compartments but we can't tell "
            "which one to use. Consider renaming your "
            "compartments please."
        )

    if most is not None:
        return most[0]
        LOGGER.warning(
            "Could not identify an external compartment by name and"
            " choosing one with the most boundary reactions. That "
            "might be complete nonsense or change suddenly. "
            "Consider renaming your compartments using "
            "`Model.compartments` to fix this."
        )
    # No info in the model, so give up
    raise RuntimeError(
        "The heuristic for discovering an external compartment "
        "relies on names and boundary reactions. Yet, there "
        "are neither compartments with recognized names nor "
        "boundary reactions in the model."
    )

cobra.Model.copartmentsが以下のどれかだと、外部環境を表す区画として認識してくれる様子。

cobra.medium.annotations.compartment_shortlist["e"] + ["e"]
['extracellular',
 'extraorganism',
 'out',
 'extracellular space',
 'extra organism',
 'extra cellular',
 'extra-organism',
 'external',
 'external medium',
 'e']

community_compartment = "out"としたら問題なく走るようになった。

結果

species = ["Bacteroides_caccae_ATCC_43185", "Lactobacillus_rhamnosus_GG_ATCC_53103"]
model_filenames = ["AGORA/CurrentVersion/AGORA_1_03/AGORA_1_03_sbml/{s}.xml".format(s=s) for s in species]

# 共培養GEMをつくる
community = create_community_model(model_filenames, species)

# Mediumの設定:Arabinogalactoseを除外した場合
medium = community.medium
for k in medium.keys(): medium[k] = 0
for ex,val in zip(dmem["Exchange reaction ID"], dmem["Uptake rate"]):
    if community.reactions.has_id(ex): medium[ex] = val
medium["EX_ala_L(e)"] = 0
medium["EX_arabinogal(e)"] = 0
community.medium = medium
op = community.optimize()
print("DMEM with Vitamin K, Hemin: {g}".format(g=op.objective_value))
print("  {s}: {g}".format(s=community.notes["species"][0]["name"], g=op.fluxes[community.notes["species"][0]["objective"]]))
print("  {s}: {g}".format(s=community.notes["species"][1]["name"], g=op.fluxes[community.notes["species"][1]["objective"]]))

# Mediumの設定:Arabinogalactoseありの場合
medium = community.medium
for k in medium.keys(): medium[k] = 0
dmem = pd.read_table("Magnusdottir-2017-Nat_Biotechnol/Medium-DMEM_6429.tsv")
for ex,val in zip(dmem["Exchange reaction ID"], dmem["Uptake rate"]):
    if community.reactions.has_id(ex): medium[ex] = val
medium["EX_ala_L(e)"] = 0
community.medium = medium
op = community.optimize()
print("DMEM with Vitamin K, Hemin, Arabinogalactan: {g}".format(g=op.objective_value))
print("  {s}: {g}".format(s=community.notes["species"][0]["name"], g=op.fluxes[community.notes["species"][0]["objective"]]))
print("  {s}: {g}".format(s=community.notes["species"][1]["name"], g=op.fluxes[community.notes["species"][1]["objective"]]))

各細菌の生育率は以下のように推定された。Lactobacillus rhamnosus GG単独のGEMにFBAを適用すると、Arabinogalactanが無いとDMEM培地では生育しないと推定されたが*3、共培養するとArabinogalactanなしでも生育すると推定された。

  • Arabinogalactose あり
    • Bacteroides_caccae_ATCC_43185: 0.4386544423746743
    • Lactobacillus_rhamnosus_GG_ATCC_53103: 0.22080916149500934
  • Arabinogalactose なし
    • Bacteroides_caccae_ATCC_43185: 0.24543405369998675
    • Lactobacillus_rhamnosus_GG_ATCC_53103: 0.11140717070316093

*1:論文ではlumenと呼んでいる

*2:https://opencobra.github.io/cobratoolbox/latest/tutorials/tutorialMicrobeMicrobeInteractions.html

*3:論文と同じ&前回記事参照

マイクロバイオーム解析④:COBRApyで細菌の代謝をモデリングする

COBRApy

以前の記事で触れたFlux balance analysis (FBA) にはMATLAB用のパッケージCOBRAが使われている論文が多い。筆者の解析環境にはMATLABはインストールされていないので*1COBRAの機能をPythonで使えるようにしたパッケージCOBRApyを使おうと思うが、COBRAと比較してチュートリアルが充実しておらず*2、実際に使う上で少しハードルが高い印象を受けた。この記事では、COBRAの開発グループが2017年にNature Biotechnologyへ掲載した論文のデータをCOBRApyを使って再現してみることで、COBRApyの機能を理解しようと思う。この論文は、ゲノムデータ・文献情報をもとに腸内細菌773種のGEMs (Genome-scale metabolic models)のセットAGORAを作成し、各細菌の代謝の特徴や、2種の細菌が共存する場合の代謝相互作用をFBAを使って推定している。

www.nature.com

本記事の検証に使ったコードはgithubにアップした。

github.com

データを準備

git clone https://github.com/tmaruy/cobrapy.git
cd cobrapy

# AGORAをダウンロード
git clone https://github.com/VirtualMetabolicHuman/AGORA.git

# メタデータをダウンロード
curl "https://static-content.springer.com/esm/art%3A10.1038%2Fnbt.3703/MediaObjects/41587_2017_BFnbt3703_MOESM6_ESM.xlsx" --output Magnusdottir-2017-Nat_Biotechnol/41587_2017_BFnbt3703_MOESM6_ESM.xlsx
curl "https://static-content.springer.com/esm/art%3A10.1038%2Fnbt.3703/MediaObjects/41587_2017_BFnbt3703_MOESM10_ESM.xlsx" --output Magnusdottir-2017-Nat_Biotechnol/41587_2017_BFnbt3703_MOESM10_ESM.xlsx

Fig. 1

Fig. 1bではAGORAに存在する各細菌のGEMsの特性 (e.g. GEMsに存在する代謝反応・代謝物の数) を示している。これを再現してみる。AGORAはGEMsをSBML形式 *3 で保存している。COBRApyではこれをcobra.Model型のインスタンスとして読み込む。cobra.Model型の基本的な特性やモデルの読み込みについてはこのノートCOBRApyのマニュアルを参照されたい。

import numpy as np
import pandas as pd
import glob
import cobra
from plotnine import * # ggplot2-like plotting in python

# GEMのリストを取得
sbml_files = np.sort(glob.glob("AGORA/CurrentVersion/AGORA_1_03/AGORA_1_03_sbml/*xml"))
bacteria_names = [f.split("/")[-1].replace(".xml", "") for f in sbml_files]

# 細菌メタデータ
s_meta = pd.read_excel("Magnusdottir-2017-Nat_Biotechnol/41587_2017_BFnbt3703_MOESM6_ESM.xlsx", skiprows=2)
s_meta = s_meta.loc[s_meta.ModelAGORA.isin(bacteria_names), :] # GEMが見つからない細菌を除外

# GEMから代謝反応、代謝物質のリストを抽出
s2r = []; s2m = []
for s in tqdm.tqdm(s_meta.ModelAGORA):
    file = "AGORA/CurrentVersion/AGORA_1_03/AGORA_1_03_sbml/{s}.xml".format(s=s)
    model = cobra.io.read_sbml_model(file)
    s2r.append(pd.DataFrame({"species": s, "reactions":[x.id for x in model.reactions]}))
    s2m.append(pd.DataFrame({"species": s, "metabolites":[x.id for x in model.metabolites]}))
s2r = pd.concat(s2r, ignore_index=True)
s2m = pd.concat(s2m, ignore_index=True)

# Fig. 1b: 代謝反応の数
figmat = s2r.groupby("species").size().rename("n_reaction").reset_index()
figmat = pd.merge(figmat, s_meta.rename(columns={"ModelAGORA":"species"}).loc[:, ["species", "Phylum"]], how="left", on="species")
(
    ggplot(figmat) +
        geom_boxplot(aes(x="Phylum", y="n_reaction")) +
        theme_bw() +
        theme(axis_text_x=element_text(angle=90), figure_size=(5,3)) +
        scale_y_continuous(limits=[500, 2000], breaks=[500, 1000, 1500, 2000]) +
        ggtitle("Number of reactions")
)

# Fig. 1b: 代謝物質の数
figmat = s2m.assign(metabolite_id = lambda x: [y[:-3] for y in x.metabolites],
                    location = lambda x: [y[-2] for y in x.metabolites]).\ # merge extra/intracellular metabolites
            groupby(["species"]).metabolite_id.nunique()
figmat = pd.merge(figmat, s_meta.rename(columns={"ModelAGORA":"species"}).loc[:, ["species", "Phylum"]], how="left", on="species")

(
    ggplot(figmat) +
        geom_boxplot(aes(x="Phylum", y="metabolite_id")) +
        theme_bw() +
        theme(axis_text_x=element_text(angle=90), figure_size=(5,3)) +
        scale_y_continuous(limits=[400, 1300], breaks=[400,600,800,1000,1200]) +
        ggtitle("Number of metabolites")
)

結果

Fig. 1bのデータと類似した傾向が認められた。

f:id:auroratummy:20220103012559p:plain f:id:auroratummy:20220103012621p:plain

Supplementary Table 6

FBAでは、培地条件など、周辺環境中の物質の組成を変えることによる影響をシミュレーションすることができる。Supplementary Table 6では、2種類の食事 (Western dietとHigh fiber diet) を摂取した場合、酸素存在化・無酸素条件下の最大生育率 *4 を細菌ごとに計算している。環境中成分の設定方法は本家のマニュアルに記載されている。有酸素条件では、このページを参考に、細菌による酸素取り込みの流束を10 mmol/gDW/h *5 とした。

ここでは、以下のようなコードでSupplementary Table 6のデータを再現してみた。全773種について計算すると時間がかかるので、ランダムに選んだ50種について、4種類の生育条件 (2種類の食事と2種類の酸素条件) での生育率を推定し、論文のデータと比較した。

growth_rate_check = {"species":[], "Western_Anaerobic":[], "Western_Aerobic":[], "HFD_Anaerobic":[], "HFD_Aerobic":[]}

for s in tqdm.tqdm(species):
    growth_rate_check["species"].append(s)
    file = "AGORA/CurrentVersion/AGORA_1_03/AGORA_1_03_sbml/{s}.xml".format(s=s)
    model = cobra.io.read_sbml_model(file)

    # Western diet
    with model:
        # Anaerobic
        medium = model.medium
        for k in medium.keys(): medium[k] = 0
        for ex,val in zip(diet["Exchange reaction ID"], diet["Western Diet"]):
            if model.reactions.has_id(ex): medium[ex] = val
        model.medium = medium
        op = model.optimize()
        growth_rate_check["Western_Anaerobic"].append(op.objective_value)

        # Aerobic
        medium["EX_o2(e)"] = 10
        model.medium = medium
        op = model.optimize()
        growth_rate_check["Western_Aerobic"].append(op.objective_value)

    # High fiber diet
    with model:
        # Anaerobic
        medium = model.medium
        for k in medium.keys(): medium[k] = 0
        for ex,val in zip(diet["Exchange reaction ID"], diet["High fiber Diet"]):
            if model.reactions.has_id(ex): medium[ex] = val
        model.medium = medium
        op = model.optimize()
        growth_rate_check["HFD_Anaerobic"].append(op.objective_value)

        # Aerobic
        medium["EX_o2(e)"] = 10
        model.medium = medium
        op = model.optimize()
        growth_rate_check["HFD_Aerobic"].append(op.objective_value)

growth_rate_check = pd.DataFrame(growth_rate_check)

結果

論文と使っているAGORAのバージョンが違うからか、最適化に使っているソルバーの違いからか、若干の差異はあるが、論文のデータとおおよそ相関する結果が再現できた。

f:id:auroratummy:20220103065051p:plain

Fig. 3

ここでは各細菌がどの物質を利用して生存できるか、どの物質を産生できるか、をFVAで解析している。FVAは目的関数がある条件を満たす場合に各代謝反応が取りうる流束の範囲を推定する方法である。詳しい解説は次の論文を参照のこと。

bmcbioinformatics.biomedcentral.com

論文によると生育率が0.001  h^{-1}を上回る状況(細菌が増加する状況)で各代謝がとりうる流束の範囲をFVAで推定しているらしい。ここではfraction_of_optimumというオプションを使って同様の条件を設定してみた。

論文と同様に、Exchange reaction *6 の最小流束が0よりも下回る場合に「細菌がその物質を環境から取り込む」、最大流束が0よりも高い場合に「細菌がその物質を環境へ分泌する」と考える。

fig3_reactions = pd.read_table("Magnusdottir-2017-Nat_Biotechnol/fig3_metabolites.tsv") # fig.3に含まれている反応のリスト

def fig3_metabolic_potential(species):
    species2potential = []
    for s in tqdm.tqdm(species):
        file = "AGORA/CurrentVersion/AGORA_1_03/AGORA_1_03_sbml/{s}.xml".format(s=s)
        model = cobra.io.read_sbml_model(file)
        op = model.optimize()
        value = op.objective_value
        fva = cobra.flux_analysis.flux_variability_analysis(model, fraction_of_optimum=0.001/value)
        fva_ex = fva.loc[[i.startswith("EX") for i in fva.index],:].\
                    reset_index().rename(columns={"index":"Reaction"}).\
                    merge(fig3_reactions, how="right", on="Reaction").\
                    loc[:, ["Reaction", "Reaction name", "Category", "minimum", "maximum"]].\
                    assign(species = s)
        species2potential.append(fva_ex)
    return pd.concat(species2potential)

結果

主要なPhylumであるActinobacteria, Bacteroidetes, Firmicutes, Proteobacteriaの代謝能を比較してみた。ここでは20種ずつ計算した。Glycerolの代謝がBacteroidetesで少なく、Succinateの代謝がFirmicutesで少ないなど、Fig3で確認できた傾向が再現できた。

f:id:auroratummy:20220103081924p:plain

またFusobacterium 16種の代謝能を調べたところ、論文で触れられているのと同様に、全種がButyrate産生のポテンシャルを有していることがFVAで推測された。

Fig. 4

Fig. 4ではBacteroides caccae ATCC 43185とLactobacillus rhamnosus GGの生育条件を調べている。ここではDMEM 6429という培地に特定の物質を加えることで生育が変わるかどうかを推定している。

dmem = pd.read_table("Magnusdottir-2017-Nat_Biotechnol/Medium-DMEM_6429.tsv")

s = "Bacteroides_caccae_ATCC_43185"
file = "AGORA/CurrentVersion/AGORA_1_03/AGORA_1_03_sbml/{s}.xml".format(s=s)
model = cobra.io.read_sbml_model(file)

medium = model.medium
for k in medium.keys(): medium[k] = 0
for ex,val in zip(dmem["Exchange reaction ID"], dmem["Uptake rate"]):
    if model.reactions.has_id(ex): medium[ex] = val
model.medium = medium
op = model.optimize()
print("DMEM with Vitamin K, Hemin, Arabinogalactan: {g}".format(g=op.objective_value))

medium["EX_arabinogal(e)"] = 0
model.medium = medium
op = model.optimize()
print("DMEM with Vitamin K, Hemin: {g}".format(g=op.objective_value))

#medium["EX_mqn7(e)"] = 0
#medium["EX_mqn8(e)"] = 0
medium["EX_pheme(e)"] = 0
model.medium = medium
op = model.optimize()
print("DMEM only: {g}".format(g=op.objective_value))

結果

論文での予測と同様にBacteroides caccae ATCC 43185はArabinogalactan存在化で生育率が上昇することが示唆された。Lactobacillus rhamnosus GGも論文と同様にNicotinic acidが生育に寄与することが示唆された。

  • Bacteroides caccae ATCC 43185
    • DMEM with Vitamin K, Hemin, Arabinogalactan: 0.5330796480552902
    • DMEM with Vitamin K, Hemin: 0.2768354439324188
    • DMEM only: 0.2768354439324188
  • Lactobacillus rhamnosus GG
    • DMEM with Vitamin K, Hemin, Arabinogalactan: 0.0
    • DMEM with Vitamin K, Hemin, Arabinogalactan, Nicotinic acid: 0.312092469574153
    • DMEM with Nicotinic acid: 0.312092469574153
    • DMEM with Nicotinic acid w/o L-alanine: 1.4612128398466525e-16

Fig. 4dで行われているような細菌間相互作用の解析も興味深いので次の記事で触れたい。

*1:そもそもMATLABを使いたくないので

*2:COBRAユースケースhttps://opencobra.github.io/cobratoolbox/latest/tutorials/index.htmlにまとめられている

*3:Systems biology markup language. xml形式. 代謝やシグナルパスウェイを記述する用途でよく用いられる

*4:細菌の増殖を目的関数とした最適化を実施

*5:concentration per gram dry weight of cells and hour

*6:細胞外と細胞内での物質交換

マイクロバイオーム解析③:bioBakery, Kraken/Bracken

bioBakery

bioBakeryはHarvardのCurtis Huttenhowerのグループが開発したマイクロバイオーム解析のツール群である。マイクロバイオーム関連の研究を見ているとよく名前が出てくるツールが多数含まれているので、改めてどのようなツールが提供されているのか、気になったものについて調べてみる。

  • MetaPhlAn: ショットガンメタゲノムデータから由来細菌種・菌叢組成を推定
  • HUMAnN: ショットガンメタゲノムデータからを機能遺伝子を同定・機能プロファイルの推定
  • PICRUSt: メタ16Sデータから菌叢の機能プロファイルを予測
  • MelonnPan: ショットガンメタゲノムデータから機械学習代謝物プロファイルを予測
  • StrainPhlAn: マーカー遺伝子中の多型情報をもとにStrainレベルの情報を取得
  • WAAFLE: メタゲノムアッセンブリで得られたコンティグから遺伝子水平伝播を検出
  • KneadData: ショットガンメタゲノムデータの前処理ツール

huttenhower.sph.harvard.edu

MetaPhlAn3

MetaPhlAnは、約17,000種(~13,500 bacterial and archaeal, ~3,500 viral, and ~110 eukaryotic)のゲノムから構築した、系統特異的な遺伝子のリストをメタゲノムデータから探索することで、細菌叢の組成を推定する方法である。細菌の全ゲノムデータへマッピングするよりもはるかに効率よく菌叢組成を推定できる。2012年に公開された最初のバージョンから徐々に遺伝子リストが更新されており、2020年にバージョン3が発表された。

www.nature.com

# データベース構築
metaphlan --install --bowtie2db ${db} # $dbにデータベースのファイルが作られる

# 細菌種同定
metaphlan \
   ${fastq1},${fastq2} \ # Input: ペアエンドの場合はカンマで区切って指定する
   --bowtie2db ${db} \ # DB
   --index latest \# DBのバージョン (default: "latest" 最新のDBを使う $dbになければインストールする)
   --input_type fastq \ # 
   
   --nproc 8 \ # CPU数
   --tax_lev {a, k, p, c, o, f, g, s} \ # 同定する系統レベル (defalt: a (all))
   --unknown_estimation \ # 未知の細菌がどの程度いるかを推定する
   --bt2_ps {very-sensitive, sensitive, very-sensitive-local, sensitive-local} \ # bowtie2のモード (default: very-sensitive)
   --min_mapq_val 5 \ # defualt: 5
   --add_viruses \ # ウイルスも検出する
   --ignore_eukaryotes \ # 真核生物を考慮しない
   --ignore_bacteria \ # バクテリアを考慮しない
   --ignore_archaea \ # アーキアを考慮しない
   --read_min_len 70 \ # この長さより短いリードは考慮しない

   -bowtie2out ${bowtie2_out} \ # 中間ファイル (bowtie2の出力) これを残しておくと再計算がすぐできる
   -s ${sam} \ # 
   -o ${result} \ #  
   -t rel_ab # return relative abundance (default: rel_ab) 

Kraken/Bracken

Kraken/BrackenはbioBakeryではないがこちらもよく目にするのでチェック。MetaPhlAnと同様にショットガンメタゲノムデータの細菌種同定・菌叢組成推定のためのツール。Krakenは系統特異的なk-merの情報を元に各リードの由来生物種をLowest Common Ancestor *1で推定する。BrackenはKrakenの後に開発された手法。Lowest Common Ancestorでは多くのリードが種・属レベルでアサインすることが困難なので、ベイズ的手法を使ってどの種にアサインされるか確率的な推定を行う。

  •  P(S_{i}|G_{j}) = \frac{P(G_{j}|S_{i})P(S_{i})}{P(G_{j})}
    •  S_{i}: Genome of Species  i
    •  G_{j}: Taxonomic group  j
    •  P(S_{i}|G_{j}):  G_{j}にKrakenでアサインされたリードが S_{i}に由来する確率
    •  P(G_{j}|S_{i}):  S_{i}由来のリードが G_{j}アサインされる確率
    •  P(G_{j}) = 1とする

 P(G_{j}|S_{i})は以下のように求める。ゲノム S_{i}からリードと同じ長さ rの配列を生成し、そのうちKrakenで G_{j}アサインされる割合を調べる。

  •  P(G_{j}|S_{i}) = \frac{N_{G_{j}(S_{i})}}{(L_{i}-r+1)}
    •  N_{G_{j}(S_{i})}:  S_{i}から生成した長さ rの配列のうち G_{j}アサインされた数
    •  L_{i}-r+1: 長さ L_{i}のゲノム S_{i}から生成可能な長さ rの配列数

 P(S_{i})は以下のように求める。ゲノム S_{i}アサインされたリードの数から推定する。

  •  U_{S_{i}} = \frac{N_{S_{i}}}{(L_{i}-r+1)}: ゲノム S_{i}から生成した長さ rの配列のうち S_{i}に直接アサインされる割合
    •  N_{S_{i}}: ゲノム S_{i}から生成した長さ rの配列のうち S_{i}アサインされた数
  •  \hat{K}_{S_{i}} = \frac{K_{S_{i}}}{U_{S_{i}}}
    •  K_{S_{i}}: 解析したメタゲノムデータで S_{i}アサインされたリードの数
    •  \hat{K}_{S_{i}}: 補正した K_{S_{i}} *2
  •  P(S_{i}) = \frac{\hat{K}_{S_{i}}}{\sum_{a=1}^{n} \hat{K}_{S_{a}}}

以上のように求めた P(S_{i}|G_{j})を用いて、 G_{j}アサインされたリードの数から、各ゲノム S_{i}アサインされるリード数を推定する。

# krakan2
kraken2 \
  --db ${database} \# k-mer database
  --confidence 0.1 \ # threshold (default: 0.0)
  --threads 8 \ #スレッド数
  --output ${kraken_output} \# Output file 1 (各リードがどの系統にアサインされたか)
  --report ${kraken_report} \# Output file 2 (各系統にLCAでアサインされたリードの数)
  --gzip-compressed \# Inputがgzipされている場合は指定する必要あり
  --memory-mapping \# k-mer databaseをメモリにロードしたくない場合
  --paired ${forward_reads} ${reverse_reads}
 
# bracken
bracken \
  -d ${database} \# k-mer db (same as kraken)
  -i ${kraken_report} \# Input from kraken
  -w ${bracken_output1} \# Output file 1 (Brackenのアルゴリズムで推定された各系統のリード数)
  -o ${bracken_output2} \# ${bracken_output1}の内容をKrakenのフォーマットで出力したもの
  -r 100 \# リード長 (上記の計算で使用する)
  -l {S, G, F} # S: species, G: genus

HUMAnN

HUMAnNは機能遺伝子の同定・菌叢の機能プロファイリングを効率的に行うパイプライン。(i) Prescreen: MetaPhlAnで細菌種特異的なマーカー遺伝子を探索 (ii) Nucleotide search: 同定された細菌種のゲノムをリファレンスとしてIndexを作成し、マッピング (iii) Translated search: マップされなかったリードをUniref90に対してDiamondでblastx-likeな検索を実施、という3ステップ。

# ペアエンドのデータを連結 (humannはペアの情報を考慮しない)
cat ${forward_fastq} ${reverse_fastq} > ${fastq}

# 
humann \
   -i ${fastq} \# Input file
   -o ${output_directory} \# 結果が出力されるディレクトリ
   --taxonomic-profile taxonomic_profile.tsv \# 
   --remove-temp-output \# 中間ファイルを消す (default: off)
   --o-log ${log_file} \# ログファイル (default: temp/sample.log)
   --output-basename ${basename} \# 出力ファイルのPrefix
   --output-format {tsv, biom} \# default: tsv

   --bypass-prescreen \# Prescreenをスキップ (Nucleotide searchには全細菌のゲノムを使う)
   --bypass-nucleotide-search \# Nucleotide searchをスキップ
   --bypass-translated-search \# Translated searchをスキップ
   
   --threads 8 \# スレッド数
   --nucleotide-database ${database}  

引用元:https://github.com/biobakery/biobakery/wiki/humann3

KneadData

ショットガンメタゲノムデータの前処理のためのツール。(i)低クオリティリードの除外と、(ii)アダプター配列の除去をTrimmomaticで、(iii)リピート領域に由来するリードの除去をTRFで、(iv)ホストゲノムに由来するリードをBowtie2で除外する。ホストゲノムのリファレンスにはhg37*3に加えて、Decoy Genomeと呼ばれるヒトゲノム由来ながらhg38にマッピングされないリードから構築されたコンティグ*4や、公共バクテリアゲノムによくコンタミしていることが報告されているヒト由来配列*5を用いる。メタトランスクリプトームデータの場合はリファレンスにはホストのトランスクリプト配列を用いる。リファレンスはコマンドkneaddata_databaseで既に構築されているものをダウンロードすることができる。

# ゲノムの場合 ($REFDIRにダウンロードされる)
kneaddata_database --download human_genome bowtie2 $REFDIR

# トランスクリプトの場合
kneaddata_database --download human_transcriptome bowtie2 $REFDIR

kneaddataは以下のコマンドで実行可能。

kneaddata \
   --input1 ${fastq_R1} \ # fastq file, forward reads
   --input2 ${fastq_R2} \ # fastq file, reverse reads
   -db ${REFDIR} \ # データベース (上述の流れで準備)

   -o ${output} \ # output
   --output-prefix ${prefix} \ # 
   --reorder \ # outputの配列をinputと同じ順番で並べる (default: off)

   -p ${process} \ # プロセス数
   -t ${threads} \ # スレッド数
   --max-memory 80g \ # 使用メモリ上限

   --trimmomatic ${path_trimmomatic} \ # Trimmomaticの場所 (default: $PATH)
   --run-trim-repetitive \ # FASTQCで出てくるoverrepresented sequencesを除外するオプション (default: off)
   --trimmomatic-options="ILLUMINACLIP:TruSeq3-SE.fa:2:30:10" \ # Trimmomaticのオプション (オプションごとに分ける必要がある)
   --trimmomatic-options="SLIDINGWINDOW:4:20" \
   --trimmomatic-options="MINLEN:50" \ 

   --bowtie2 ${path_bowtie2} \ # Botwie2の場所 (default: $PATH)
   --bowtie2-options ${bowtie2_option} \ # default: --very-sensitive-local
   --decoantaminate-pairs {strict, lanient, unpaired}\ # strict: forward/reverseの片方でもマップされれば除外 (default)

   --run-trf # tandem repeat に対応する配列を除外 (default: on)

引用元:https://github.com/biobakery/kneaddata

*1:あるリードが細菌Aと細菌Bに対応するk-merを含有する場合はAとBに共通する系統をアサインする、例えばAとBが同じ属であれば、その属をアサインする

*2:S_i由来のリードのうち1割しかS_iにアサインされないのであれば、実際のS_i由来のリード数はK_S_i / 0.1と推定できる

*3:hg38ではない

*4:詳細はhttp://www.cureffi.org/2013/02/01/the-decoy-genome/、ヒトゲノムに組み込まれたウイルス由来の配列などが含まれる

*5:https://pubmed.ncbi.nlm.nih.gov/31064768/