ショートリードからの高分解能菌株レベルマイクロバイオーム組成解析

2023年8月18日 11:20

オープンアクセス
公開日：2023年8月17日
ショートリードからの高分解能菌株レベルマイクロバイオーム組成解析

https://microbiomejournal.biomedcentral.com/articles/10.1186/s40168-023-01615-w

ヘルイ・リャオ、ジー・ヨンシン、ヤンニ・スン
マイクロバイオーム第11巻、記事番号：183（2023）この記事を引用する

4 Altmetric

メトリクス詳細

要旨
背景
微生物群集の動態を理解する上で、菌株レベルの組成解析は重要なステップである。メタゲノムシーケンスは、宿主関連サンプルや環境サンプル中の微生物組成をプローブするための主要な手段となっている。組成解析ツールは数多く存在するが、菌株レベルの解析における課題、すなわち、類似性の高い菌株ゲノムや、1つのサンプル中に1つの生物種に属する複数の菌株が存在することに対処するために最適化されていない。そこで本研究では、ショートリードを用いた高分解能かつ高精度な菌株レベル解析ツールを提供することを目的とした。

結果
本研究では、菌株同定精度と計算の複雑さのバランスを取るために、新しいツリーベースのk-mersインデックス構造を採用したStrainScanと名付けた新しい菌株レベル構成解析ツールを発表する。我々は、StrainScanを多数のシミュレーションデータと実際のシーケンスデータで広範囲にテストし、Krakenuniq、StrainSeeker、Pathoscope2、Sigma、StrainGE、StrainEstなどの一般的な菌株レベル解析ツールとのベンチマークを行った。その結果、StrainScanはひずみレベルの組成解析において、最先端のツールよりも高い精度と分解能を持つことが示された。StrainScanは、菌株レベルでの複数菌株の同定において、F1スコアを20％向上させた。

結論
新規のk-merインデックス構造を使用することにより、StrainScanは既存のツールよりも高い分解能で株レベル解析を提供することができ、1つのサンプルまたは複数のサンプルにおいて、より有益な株組成解析を返すことができる。StrainScanは、ショートリードと参照株のセットを入力とし、そのソースコードはhttps://github.com/liaoherui/StrainScan で自由に利用できる。

ビデオ

背景
ゲノム変異により、1つの生物種内の菌株が異なる代謝および機能的多様性を持ちうることを示す証拠が蓄積されつつある[1,2,3]。同じ生物種に属する菌株でも、高い配列多様性と異なる遺伝子構造を示すことがある[4]。菌株にユニークな遺伝子やSNPがあれば、新たな酵素機能、抗生物質耐性、病原性、異なる感染ウイルスなどをもたらす可能性がある。例えば、大腸菌には少なくとも数千の菌株が同定されており、その中には病原性因子を持つものもあれば、常在菌のものもある。顕著な例として、2011年にドイツで発生したO104:H4株による大腸菌アウトブレイクが挙げられる。O104:H4株は志賀毒素をコードするプロファージとその他の病原性因子を獲得していた[5]。

菌株が異なれば生物学的特性も異なるため、菌株を特定することはマイクロバイオームの構成と機能解析の両方にとって重要である。メタゲノム解析データは、宿主に関連したサンプルや環境サンプルから得られた遺伝子の配列決定データを含んでおり、細菌の菌株レベルの組成を研究するための主要な情報源となっている。さまざまなサンプルにおける菌株の遺伝子型や表現型に関する新たな知識を生み出す研究が増えている。例えば、Pollardらは、198の海洋メタゲノムにおいて、多くの一般的な細菌種が地理的な場所に関連した菌株レベルの組成を持つことを示した[6]。密接に関連した研究では、クローン病患者の腸内マイクロバイオームにおいて、優勢な大腸菌株が経時的に変化することが示された[7]。ヒトの腸内にごく普通に存在するもう一つの細菌であるP. copriは、その菌株と宿主の地理的位置や食習慣との間に密接な関係があることが証明されている [8, 9]。潜在的なプロバイオティクスであるA.ムチニフィラの菌株の中には、抗炎症作用を持つものがあり、肥満や糖尿病に有益な効果をもたらす可能性があることが分かっている [10]。さらに、人体の各部位における菌株の分布には違いがある。例えば、過去の研究[11]では、身体の異なる部位から採取されたC. acnesとS. epidermidisの菌株は、不均一で多系統であることが判明している。

図1
図1
StrainScanは、同定されたクラスター内の菌株を検索することで、より高い菌株レベルの分解能を実現する。対照的に、StrainGEやStrainEstのようなクラスター・レベルのツールは、同定されたクラスターの代表菌株のみを返し、クラスター内の他の菌株は検索しない。「S1」と「S2」は2つの入力メタゲノミックサンプルである。B-C 216株の大腸菌を含む実際のクラスター（クラスター1と命名）の菌株間の遺伝子含量の違い。Bの10株は合計1722個の株特異的遺伝子を持つ。Cの "GCF_001695515 "と "GCF_013167975 "は、それぞれクラスター1で最も長い株と最も短い株である。

フルサイズ画像
株レベルの解析の重要性にもかかわらず、種レベル以下の分類学的解析を行うことは依然として困難である。その1つの課題は、1つのサンプル中に類似性の高い複数の株が同時に存在しうるという事実に由来する[12]。例えば、最近のある研究 [13] では、ヒトの糞便サンプル中に2～3株の表皮ブドウ球菌が共存し、そのMash [14] 距離は約0.005であることが判明した。同様に、ヒトの皮膚マイクロバイオームの重要な構成要素であるC. acnesの複数の菌株が、しばしば複雑な混合物を形成することを示す報告もある [15]。これらの共存菌株の中には高い配列類似性を示すものがあり、マッシュ距離は約0.0004である。さらに、2144のヒトの糞便メタゲノムを解析した研究[16]では、多数のサンプルに非常に類似したバクテロイデス・ドレイ株が共存していることが明らかになった。一般的に使用されているメタゲノム解析ツールは、異なる菌株を区別するようには設計されていない。菌株解析ツールはあるが、同じ集団から複数のサンプルを必要としたり[17]、優占株のみを出力したり[18,19,20,21]、菌株間の類似性に制約を与えたりする[22]。すぐに続く第2の課題は、菌株レベルの同定の分解能である。ここでの分解能は参照データベースのサイズに反映され、参照菌株の数が多いほど分解能が高いことを示す[23]。いくつかの菌株は非常に高い類似性を共有しているが、遺伝的差異を無視できる類似性のカットオフ値は知られていない。例えば、病原性大腸菌CFT073とプロバイオティクス大腸菌Nissle 1917の配列類似度はそれぞれ99.98%である[24]。同様に、あるファージ-宿主共進化研究[25]では、細菌株が高いゲノムANI（99.9％以上）を含んでいても、菌株が異なるファージに感染し、異なる防御機構や吸着機構を示す可能性があることがわかった。菌株レベルの多様性が高い一部の菌種では、わずかなSNVが表現型の変異につながる可能性がある[26, 27]。従って、より高い解像度で遺伝子型と表現型の関係をより正確に特徴付けることができる。StrainGE[28]やStrainEst[29]などのツールは、菌株の混合を解きほぐすように設計されていますが、サンプリングされた菌株ゲノムデータベースの代表的な菌株を報告することに限定されています。これらのクラスタリングのカットオフ値（0.9 k-mer Jaccard類似度（StrainGE）または99.4% ANI（StrainEst））でも、細菌によっては大きなクラスタになることがあります。StrainGEは、サンプル中の同定された代表的な菌株に対するSNP/欠失をさらに同定することができますが、同定されたクラスター中の特定の菌株を特定することはできません。2つのk-merベースのツール、Krakenuniq [30]とStrainSeeker [31]も、データベース内の菌株が高い類似性を共有している場合、菌株レベルの同定では分解能が非常に低い。第三の課題は、低存在株の同定である。例えば、de novo株構築ツール[32, 33]は、アセンブリベースの戦略を用いて株を再構築することを目的としているが、正確な株の再構築を達成するためには、通常、高いカバレッジの株が必要となる。さらに、多くの菌株解析ツール[34,35,36]も、正確な同定を行うためには10倍以上の菌株カバー率を必要とする。そのため、これらのツールでは、カバレッジの低い菌株を同定することが課題となっている。最後の課題は菌株同定に要する時間である。最近発表された研究[23, 37]によると、Sigma[38]やPathoscope2[39]などのアライメントベースの菌株レベル同定ツールのほとんどは、データベースが大規模な場合、計算コストが高くなる可能性がある。大規模な参照データベースは、種内多様性のカバレッジを高めることができる一方で、より多くの計算資源を必要とする。

したがって、メタゲノムデータに対して、より高感度、高精度、効率的な株レベル解析を提供することが急務である。本研究では、メタゲノムデータや全ゲノムシーケンスデータを含むシーケンスデータから既知の菌株を正確に検出できるオープンソースツールStrainScanを紹介する。解像度と計算の複雑さのバランスを取るために、我々は、通常、不均一な類似性分布を示す多数の菌株に対して、新しい階層的k-mersインデックス構造を開発した。まず、類似度の高い菌株をクラスターに分類する。次に、クラスタ検索のためのツリーベースのインデックス構造である新規なクラスタサーチツリー（CST）を設計する。各ノードに含まれるk-merの数を注意深くバランスさせることにより、CSTを最適化し、存在量の少ない菌株の誤同定を防ぐ。第2ステップでは、菌株固有のk-merとSNVや構造変異を表すk-merを用いて、どの菌株が存在する可能性が高いかを判定する。StrainScanの最終出力には、同定された菌株とその存在量が含まれる。同定されたクラスター内の菌株を検索することにより、StrainScanはStrainGEやStrainEstのようなクラスター・レベルのツールよりも高い解像度を達成する。図1Aに示すように、解像度が異なると、観察結果や結論も異なってくる。StrainScanはサンプルS1に含まれる2つの異なる菌株を同定できるが、StrainGEやStrainEstでは、同じクラスターの菌株であるため区別できない。同様に、2つのサンプル（S1とS2）を比較する際、クラスターではなく特定の菌株をピンポイントで特定することで、より正確な遺伝子組成に基づく解析が可能になる。

StrainScanと他の利用可能なツールとのベンチマークを複数のシミュレーションおよび実際のシーケンスデータセットで行うことで、StrainScanが最先端のツールよりも高い精度で株レベルの組成を出力できることを実証した。特に、StrainGEのような最先端のツールと比較した場合、StrainScanは菌株レベルでの複数菌株の同定において、F1スコアを20％以上向上させた。StrainScanは、対象とする菌株のリファレンスゲノムを提供する必要がある、標的株構成解析ツールである。任意の参照ゲノムセットに対してインデキシング構造のカスタマイズ構築をサポートすることで、あらゆる細菌に適用することができる。

方法
StrainScanの概要
StrainScanは、ショートリードから既知の菌株を直接同定するように設計されています。メタゲノムデータには多くの種レベルの組成解析ツールがあるため、StrainScanへの入力は、"fastq "フォーマットのショートリードと "fasta "フォーマットの対象細菌の菌株ゲノムである。菌株同定の分解能と計算コストのバランスを取るために、高速だが粒度の粗いクラスターサーチツリー（CST）と、低速だが粒度の細かいクラスター内の菌株同定ストラテジーを組み合わせた階層的インデックス作成法を設計した。図2のフローチャートに示すように、まずクラスタツリーベースのインデックス構造を作成します。このツリー上での効率的で正確なクラスタ探索法により、まずサンプルに存在するクラスタをピンポイントで特定することができます。次に、慎重に選択したk-merを使用して、同定されたクラスター内の異なる株を区別する。階層的手法にはいくつかの利点がある。第一に、菌株間の不均一な類似性分布に対応でき、ある菌株は他の菌株よりはるかに高い類似性を共有する。類似度の低い菌株は、高速CST検索戦略によって迅速に同定できる。そして、類似度の高い菌株だけを第2ステップでより細かく区別する必要がある。第二に、階層的手法は、より多くのユニークなk-mersを使用できるようにすることで、検索精度を向上させることができる（補足表S1）。クラスタ間で共有されているk-mersは、クラスタ内検索に利用できる。第三に、階層的手法はメモリフットプリントを削減できる。階層化法を用いない場合、多数のk-mersを含むすべての文献から株を検索する必要がある。CST検索によって特定されたクラスターがあれば、StrainScanは特定されたクラスターの中で、より少ない系統とk-merを含む系統だけを検索すればよい。例えば、クラスタリング前の大腸菌参照セットのk-mers総数は192,325,016であるのに対し、クラスタリング後の最大クラスターのk-mers数は16,071,080に減少する（補足表S1）。

図2
図2
StrainScanの概要。(a) 菌株ゲノムのクラスタリングプロセスのスケッチ。対象となる細菌の株ゲノム（G1、G2、...）が与えられると、Dashing [40]を用いて全対全kマースJaccard類似度が計算される。その後、ゲノムは単結合階層クラスタリングを用いてクラスタリングされる。デフォルトでは、クラスタリングの閾値はJaccard類似度0.95に設定されている。この例では、赤の破線で表されるカットオフが与えられ、クラスタリングプロセスによってC1からC5までの5つのクラスタが出力される。(b) クラスタが与えられたら、後のクラスタレベル識別のために階層クラスタ木を構築する。(c) 同一クラスター内で異なる系統を区別するのに役立つk-merを抽出するために、コリニアブロックを生成する。(d) ステップdで、リファレンスゲノムのインデキシング構造プロセスを終了する。(e)および(f) インデキシング構造とシーケンスデータ（リード）を入力し、株を検索する。(e) クラスターを検索する。(f)反復行列乗算によって株が同定され、最終的にelastic net regressionによって相対存在量プロファイルが推定される。

フルサイズ画像
図3
図3
A CSTベースのインデックス構造におけるk-mersの割り当ての例。各ノードはそのルートとなるサブツリーに固有のk-mersを持ち、サブツリー内のほとんどの系統で共有されている。この例では、特定の色の棒がk-mersを表し、各ノードに1つずつ固有のk-mersが割り当てられている。B ノードvのk-mers集合を構築するとき、すべての葉ノードはPRE_v、SUB_v、EXT_vという3つのグループに分けられる。

フルサイズ画像
クラスターサーチツリー（CST）の構築
同一種の多数の菌株ゲノムが与えられた場合、まずアラインメントを必要としないk-mersベースの手法Dashing [40] (k = 31)を用いてJaccard類似度行列を計算した。次に、このマトリックスに基づいて凝集型階層クラスタリング（単一結合）を行い、菌株をデンドログラムにグループ化した。最後に、固定の高さカットオフH（デフォルトでは0.95）を選択し、デンドログラムを多数のクラスターに分割した。各クラスタ内の菌株はk-merベースのJaccard類似度( \ge 0.95)を有し、これは99.89%[28]の平均ヌクレオチド同一性(ANI)にほぼ相当する。

菌株が含まれるクラスターをピンポイントで特定するために、クラスターとデンドログラムをCSTに変換し、正確で効率的なクラスター検索の両方をサポートします。CSTは、各クラスタがツリーのリーフ・ノードで表される以外は、デンドログラムと同じツリー・トポロジーを維持する。さらに、各ノードとその親（または子）間の距離が、それらのJaccard類似度に関係なく一様になるように、デンドログラムの距離情報を破棄します。したがって、CSTは完全なバイナリ・ツリーとなる。クラスタ探索をサポートするために、各ノードは、このノードがルートとするサブツリーに固有なk-mersの集合を含む。k-mersマッチを行うことで、CSTは、1つまたは複数のリーフ・ノード（すなわち、クラスター）に到達するまで、左の子または右の子のいずれかを取るようにガイドする。まず、各ノードにどのようにk-mersを割り当てるかを説明する。

ノードのkマー割り当て
CSTは、木のトポロジーと各ノードに割り当てられたk-mersセットという2つの要素によって定義される。このセクションでは、クラスタ探索をサポートするために、ノードにどのようにk-mersを割り当てるかを説明する。CSTのノードvに対して、vをルートとする部分木をT_vとする。vのkマー割り当てには2つの基準がある。第一に、k-mersはT_vのリーフノードのほとんどの系統で共有されるべきである。第二に、k-mersはT_vの系統に固有である。この2つの基準を、図3Aの例を使って視覚化した。

2つの基準に従って、まずリーフノードに、対応するクラスター内の菌株から抽出したk-merを割り当てる。複数の系統が存在するクラスターでは、その系統で比較的よく保存されている特徴を表すk-merを使用するため、少なくとも˶‾alpha‾%に出現するk-merのみを残す。また、"alpha "が大きい場合は、多くの系統で共有されているk-merのみをCSTの構築に使用することを示し、"alpha "が小さい場合は、系統固有のk-merを使用することを示す。実験では、様々な αを使用してクラスター同定性能を比較しました。補足図S1の経験的な結果に従い、我々はデフォルトの㊟αを90に設定した。

以下、リーフノードvの初期kマース集合を「 \mathbf {L_v}」とする。次に、葉ノードから始めて、2つの兄弟ノード間の共有k-mersを再帰的に親ノードに向かって移動する。最後のステップでは、2つ以上のノードに出現するすべてのkマーを削除する。このプロセスの最後に、各ノードv（内部ノードまたはリーフノード）には一意なkマーの集合が含まれる。具体的には、ノードvのK_vは、式(1)に示すように集合演算を用いて構築することができる。深さd_vを有するノードvに対して、図3Bに示され、以下に定義されるように、全てのリーフノードは、vとの関係に基づいて3つのグループに分割される。

図
\begin{aligned} K_{v} = ┣L_{i} - ┣L_{i} - ┣L_{i L_{i} - ￤ビッグキャップ￤リミット {iin EXT_v} L{i} - ￤ビッグキャップ￤リミット {iin PRE_v} L{i｝ \end{aligned}
(1)
ここまでに構築されたCSTはStrainSeeker[31]で構築されたツリーに似ている。ユニークなk-mers ↪Mathbf {K_v}を使用することで、菌株クラスターを特定するための探索をガイドすることができますが、重大な限界は、一部のノードが少数のユニークなk-mersしか含まないことであり、これは多くのk-mersを含むノードよりも偽陽性（FP）マッチの可能性が高くなります。これはStrainSeekerを適用したときに観察されました。P. copri 112株から構築されたStrainSeekerデータベースを例にとろう。222個のノードのうち、21個のノードは空であり、104個のノードはk-mersが1000未満である。k-mersセットが小さいノードは偶然マッチングされる傾向があり、その結果FP同定につながる。この限界に対処するため、クラスタ探索に曖昧さを加えないk-mersを追加することで、これらのノードを増強する。CST最適化手法の詳細は補遺1.1節を参照されたい。

CSTにおけるクラスター探索
入力シーケンスデータが与えられたら、まずCSTから全てのk-mersを抽出し、Jellyfish [41]を用いて全てのショートリードに対して高速k-mersマッチを行う。次に、すべてのリードのk-mersマッチカウントをCSTにマッピングする。各ノードvには、1次元数値ベクトル ↪Mathbf {C_v} = (c_1, c_2, ..., c_{|Kv|})が割り当てられ、各セルにはk-mersマッチカウントが記録される。クラスタ探索アルゴリズムはBreadth-first Search (BFS)に基づいており、ルートから開始し、ノードのk-mersマッチをレベルごとに調べる（図4）。各ノードvのk-mersマッチ・ベクトルC_vは、二項検定に基づいてvの子孫を探索するかどうかを決定するために使用される。最終的な検索結果には、シーケンスデータに存在する株クラスターを表す1つまたは複数のリーフノードが含まれる。

スコアリングメトリクス
検索がノード v を訪問するとき、どの子ノードを訪問するかを決定するために、2つのスコアリング指標が計算される。図4に示すように、1つ目の指標は、K_v中のk-mersがシーケンシングデータに存在する割合を表す、マッチするk-mersの割合（ \mathbf {frac_v}）である。これは次のように定義される：

\begin{aligned} frac_v = \frac{|C^+_v|}{|C_v|} 。\end{aligned}
(2)
ここで、C^+_vはC_vからの全ての正のk-mersカウントを含むベクトルを表す。

2つ目の指標は平均k-mersマッチングカウント( \mathbf {abund_v})であり、これは正のマッチングカウントを持つk-mersのみを用いて計算される。また、frac_v < 0.1のとき、abund_vは0に設定される。

\また、frac_v < 0.1のときabund_vは0となる。\を0に設定します。\ 0, &{} {frac_v < 0.1｝ \ 0, &{} {frac_v < 0.1 \right。\end{aligned}
(3)
図4
図4
クラスタ探索プロセスの例。つのスコアリング指標（abund_i, frac_i）の値が各ノードvの横に示されている。検索結果にはクラスター4と5が含まれ、推定存在量はそれぞれ23.8と9である。ノード6と3は二項検定に合格しなかったので、それらの子孫をトラバースすることができなかった。Path_4とPath_5のノードはそれぞれ青とオレンジで色分けされている。K_4はL_5と同じオレンジ色のk-mersを共有している。したがって、正確なスコアリング指標を計算するために、クラスター5の推定アバンダンスA_5に基づいてC_4を調整する必要がある。

フルサイズ画像
探索戦略
vの2つのスコアリング指標を計算した後、CST[18]の探索順序を決めるために二項検定を行う。シーケンシングエラーはk-mersマッチを引き起こす可能性があるため、二項検定の主な目的は、シーケンシングエラーによるランダムマッチと、真の株の真のマッチを区別することである。ノードvが与えられたとき、abund_vが配列決定エラーによって生成されたという帰無仮説を棄却できるかどうかを調べます。

具体的には、まずabund_vとabund_p（pはvの親ノード）を最も近い整数abund_v'とabund_p'に丸めます。次に、配列決定エラー率e（デフォルトでは1％）が与えられたとき、abund_vが配列決定エラーから生成される確率が˶β（デフォルトでは0.05）より小さいとき、エラーが原因の帰無仮説を棄却する。この確率は

\begin{aligned}で推定される。P_{Xtilde{B}(abund_p', 1-e)}(Xle abund_p'-abund_v') ￤end{aligned}で推定する。
(4)
ここで、B(abund_p, 1-e)はabund_p''の試行と成功率1-eの二項分布の確率質量関数である。帰無仮説が棄却されないということは、低カバレッジのk-mersマッチとシーケンスノイズを区別できないことを示している。従って、abund_vはシーケンシングエラーによるものだと考え、vの子孫の探索を停止する。そうでなければ、帰無仮説の棄却に成功した場合、T_vに含まれる1つまたは複数の株クラスターがシーケンシングデータに存在すると考える。そして、CST探索は、vの2つの子ノードをBFSキューの末尾に追加し、後でそれらをトラバースする準備をする。従来のバイナリ探索木（BST）とは異なり、同じ親の2つの兄弟ノードは、エラーが原因の帰無仮説を両方とも棄却することができます。したがって、それらの子孫をすべてトラバースすることができ、CSTの探索結果はおそらく複数になる。

クラスターの識別
リーフノードに到達したら、リーフノードと、リーフノードから移動したk-mersを含むその祖先ノードの両方を用いて、k-mersマッチングとアバンダンス推定統計量をさらに調べる。リーフノードが1つしか識別されない場合、ルートからvまでのパスに沿ったすべてのノードを使用してk-mers統計量を計算することができる。しかし，識別されたリーフノードが複数ある場合は，すべての祖先ノードを使用すべきではない．その代わりに、リーフ・ノードvに一意に寄与するものだけを最終的なアバンダンスを計算するために使用すべきである。これらのノードは集合的にパスPath_vを構成することができ、Path_v上のすべてのk-mersマッチはリーフノードvの株にのみ由来する。Path_vを同定するために、まず、同定されたノードとしてvのみを含む最大部分木を同定する。そして、Path_vはこのサブツリーのルートからvまでのパスに相当する。クラスターの存在量を推定するために、Path_v上のノードのk-mersカウントを一括して使用すると、単一のリーフノードを使用するよりも高い信頼性が得られる。図4のPath_4を例にとると、2つのリーフ・ノード4と5が識別される。この場合、Path_4はノード8のルート付き部分木T_{v_8}のvを含むルートからリーフへのパスである。続いて、Path_v上の全てのk-mers一致数C_iを集めて、加重平均一致k-mers割合(weighted average fraction of matching k-mers) \mathbf {F_v}と、加重平均k-mers一致数(weighted average k-mers match count)¤A_v} を計算することができる：

\begin{aligned}とする。F_v = \frac{sum |C^+_i|}{sum \limits {iin Path_v} |C_i｝ \end{aligned}
(5)
\begin{aligned} A{v} = \frac{sum \limits _{i} in Path_v} |C^+_i|} (5){begin{aligned}。|C^+i|cdot abund{i}}}{sum \limits _{iin Path_v}. |C^+_i|} \end{aligned}
(6)
F_vが所定のカットオフ値（デフォルト値は0.4だが、異なる条件に適応するためにユーザーが値を変更できる）より大きい場合、vのクラスターがシーケンスデータに存在するとみなす。CST検索が終了すると、同定されたすべてのクラスターが推定アバンダンス（A_vで計算）とともに出力されます。また、シークエンシングデータが異なるクラスターに複数の株を含んでいる場合、弱いノードのオーグメンテーションの際にk-merが追加されるため、いくつかのFPが導入される。この問題に対処するための詳細な方法は、補足セクション1.2に記載されている。

クラスター内のひずみ識別
CSTは、所定のカットオフ値以下の類似性を持つクラスターを区別するために最適化されている。類似性の高い系統を区別するためにCSTを使用すると、共有k-merの割合が多いため、増強できない弱いノードが大量に発生する可能性がある。したがって、いったんクラスターを特定したら、類似度の高い系統を区別するためのきめ細かい方法が必要になる。一旦クラスターを特定すれば、区別すべき系統の数は元の問題空間に比べて大幅に減少する。従って、識別力のあるすべてのk-merを使用する余裕がある。最初の特徴量は、株特異的な領域に由来するユニークなk-merであり、ここではこれを株特異的k-merと呼ぶ。2つ目の特徴はグループ特異的k-merで、これはいくつかの株に共通する構造変異（SV）に由来する可能性がある。最近の研究[42]では、異なる株を区別するためにSVが用いられている。この研究にヒントを得て、いくつかの株に共通するSVからグループ特異的k-merを抽出する。しかし、株特異的k-merやグループ特異的k-merだけに頼っていると、場合によってはまだ解像度が低いという問題がある。例えば、図5ではStrain4とStrain5が同じグループ特異的k-mersを持っており、Strain5の株特異的k-mersがサンプル中に存在しない場合、2つの株を細かく区別することができない。そこで、さらに分解能を向上させるために、全ゲノムに存在するコアゲノム領域のSNVやインデルを含むジョイントk-mersセットを追加した[29, 38, 43]。図5に示すように、すべてのジョイントk-mersについて、各k-mersは系統特異的ではないが、各系統のジョイントk-mersセットはユニークである。しかし、ジョイントk-mersの数は最初の2種類のk-mersほど多くないことが多い（補足表S2）。同定の分解能を向上させるためには、これらを組み合わせる必要がある。これら3種類のk-merを利用することで、同定の分解能を向上させ、同時に探索空間を縮小することができる。

図5
図5
株特異的k-mers、グループ特異的k-mers、結合k-mersを用いて、同一クラスター内の5株を識別する。各菌株は固有のk-mersの組み合わせを持っている。

フルサイズ画像
これらのk-mersを効率的に抽出するために、Sibeliaz[44]を利用します。Sibeliazは、近縁ゲノム中の局所的に共線的なブロックを同定するために設計された効率的なツールです。Sibeliazによって生成されたブロックを基に、ハッシュベースのアルゴリズムを開発し、株ゲノムからこれらのk-mersを抽出し、後で使用するために行列に保存します。アルゴリズムの主な擬似コードを補足1.3節に示す。アルゴリズムへの入力は、同じクラスター内の菌株ゲノムと、Sibeliazによって生成されたブロックです。効率的なハッシュテーブルを使用することで、アルゴリズムは目的のk-merを迅速に抽出することができる。最後に、クラスターから抽出されたすべてのk-mersは、サイズM×N倍の行列Xに保存される。ここで、Mはk-mersの数、Nはこの行列内の菌株数である。複数のクラスターが存在する場合は、それぞれ対応する複数の行列が作成される。

選択されたk-merを用いた株の同定
前のステップでk-merを抽出した後、これらの特徴を菌株識別に使用する必要がある。同じクラスター内の近縁の菌株からなる複雑なコミュニティを分離するために、反復行列乗算を適用して共存するすべての菌株を決定し、弾性ネット回帰を用いてそれらの相対的存在量を予測する。

反復行列乗算の主な目的は、k-mers行列の3種類のk-mers（図5）を用いて同一クラスター内の菌株を決定することである。この目標を達成するために、QuantTB [37]と同様の反復戦略を用いて、サンプル中のk-mersとk-mers行列X中のk-mersを比較する。その方法は以下の通りである。木探索によって選択されたクラスターと、そのk-mers行列Xからのk-mersが与えられたら、Jellyfish [41]を適用して、配列決定データ中のこれらの選択されたk-mersをすべてカウントする。Jellyfishから選択された全てのk-merの出現をベクトルyとする： y= (y_{1}, y_{2}, y_{3}, ..., y_{M})^{textrm{T}}、ここでy_{i}は0であり、出現を表す。\ここで、y_{i}は0であり、行列中のi番目のk-mersの出現を表す。しかし、他のクラスターと重複するk-mersは、誤ったk-mersのマッチングや間違ったアバンダンスの推定につながる可能性がある。他のクラスターの影響を除去するために、もしあるk-mersがツリー検索によって検出された他のクラスターで見つかった場合、その出現は0に置き換えられる。 M倍N個のk-mers行列Xに対して、そのj番目の列X[ : , j]は次のように定義される：

\begin{aligned}と定義される。X[:,j] = (X[1,j], X[2,j], X[3,j], ..., X[M,j])^{{textrm{T}}, ⑬j = 1, ..., N ⑭end{aligned} と定義する。
(7)
Xとyに基づき、反復行列乗算を用いることで、試料中の可能性のある全ての系統を正確かつ迅速に検出することができる。Xとyが与えられると、この関数は各菌株についてスコアf_j= X[:,j] ㎟cdot yを計算する。なお、5パーセンタイルと95パーセンタイルを超える値は外れ値とみなし、外れ値はすべて0とします。ランク付けの後、関数はランク付けされたリストの上位1株を出力し、次に同定された株に含まれるすべてのk-merの出現を0に置き換えてyを更新する。デフォルト値は31*40=1240k-mersで、0以外の値を持つk-mersの出現数が与えられた閾値を下回るまで、各反復でスコアを計算し、最も可能性の高い株を特定し続ける。本研究ではすべてデフォルトのカットオフ値を用いて実験を行った。

サンプルに含まれる可能性のある菌株を把握した上で、エラスティックネット回帰モデルを用いて、同定された菌株の配列決定深度と相対存在量を予測する。Lassoモデルは菌株数を過小評価する傾向があり、再現性の低下につながるため、Lassoモデルの代わりに弾性ネットモデルを選択します。反復行列乗算の後、フィルタリングされたk-mers行列X'=M ゙times N'を得る。回帰係数である配列深度は、エラスティックネットのペナルティ付き残差平方和を最小化することにより予測される：

\begin{aligned}とする。\を最小化することにより予測する。(8){left}{y-X'}{N'}{beta_k>=0}{textrm{argmin}}+lambda({alpha}{left}{beta_right}{1+frac{(1-alpha )}{2}left}{beta_right}^2}){end{aligned}となる。
(8)
\αとλはモデル性能に影響する重要なパラメータであり、チューニングが必要である。クロスバリデーションに基づき、最も予測誤差の少ないモデルを得るために、"alpha "と "lambda "を調整する関数を設計しました。この最良のモデルが与えられたとき、モデルの回帰係数'alpha'を正規化することにより、系統相対存在量a=(a_1, a_2, a_3, ..., a{N'})を計算する。しかし，木探索によって複数のクラスターが検出された場合，1つの系統iの相対存在量はクラスターの存在量に応じて再計算される．そのため、各菌株iの最終的な相対存在量(RA)は以下のように計算される：

\begin{aligned}とする。RA_i = \frac{a_iC_i}{sum _{j=1}^na_jC_j}。\end{aligned}
(9)
ここで、Cは木探索によって予測されたクラスターの存在量、nは同定された全菌株の総数である。

表1 全実験の概要。グレーのブロック：シミュレーションデータ、青のブロック：モックデータセットまたはスパイクデータセット、オレンジのブロック：実際のシーケンスデータセット。すべてのデータセットの菌株構成は、シミュレーションプロセスまたは原著論文[9,15,28,45,46,47,48,49,50,51,52,53,54]によって提供されている。
フルサイズの表
予測精度の評価
各手法の性能を検証するため、テスト・カテゴリーごとに、リコール、精度、F1スコアを計算した。真陽性（TP）とは、正しく同定された菌株の数を示す。偽陰性（FN：False Negative）とは、検体中に存在するがツールによって検出されなかった菌株の数を指す。偽陽性（FP）とは、誤識別された菌株の数。

\begin{aligned}（開始精度 = ╱frac{TP}{TP+FP},╱リコール = ╱frac{TP}{TP+FN}. \end{aligned}。
\begin{aligned}とする。\ F1 = \frac{2PrecisionRecall}{Precision+Recall}. \end{aligned}
すべての実験において、Jensen-Shannon divergence（JSD）[55]を用いて、真の相対存在量と予測相対存在量の間の距離を測定した。予測された存在量と真の存在量が異なる次元を持つ場合、次元の低い方にゼロを加えてJSDを計算する。2つの確率分布TとPがあるとすると、そのJensen-Shannon divergenceは[0, 1]の間の値であり、次のように定義される：

\と定義される。JSD(T) = ￢D(T) + ￢D(P) ￢end=aligned
ここで

\begin{aligned} K = \frac{1}{2}(T + Γ P) Γend{aligned} とする。
であり、D(T | | K)はTからKへのKullback-Leibler divergenceと呼ばれ、次式で定義される：

\と定義される。D(T ||K) = ∕T(i) ∕log {frac{T(i)}{K(i)}}. \end{aligned}
結果
StrainScanは既知の菌株を同定することに重点を置いているため、菌株レベルの解析に計算上の難題をもたらす可能性のある6つの細菌についてStrainScanの性能をテストした。選ばれた細菌はすべて、少なくとも100株の配列が決定されている。その中には、大腸菌や表皮菌のように多数の既知菌株を持つものもある。結核菌のように、配列の類似性が極めて高い菌株もある。さらに、A. muciniphila、P. copri、C. acnesなど、通常ヒトの腸や皮膚など異なる生態系に生息する細菌を選んだ。StrainScanを評価するために複数の実験を行った。すべての実験の概要を表1にまとめた。まず、シミュレーションデータとスパイクされたメタゲノムデータにおいて、StrainScanが1つの菌株と複数の共存菌株を同定する能力をテストした。株の類似度や株のシーケンス深さなどのパラメータを設定して異なるデータセットを作成し、困難なシナリオにおける異なるツールの性能を比較するのに役立てた。次に、StrainScanを3つの模擬コミュニティデータセットでテストしました。これにより、菌株構成が既知の実際のシーケンスデータで異なるツールを評価することができます。第三に、様々な深さの94の実シーケンスデータセットでStrainScanをテストした（補足表S3）[21, 29, 37, 39, 52]。実際のシーケンスデータには通常、菌株構成のグランドトゥルースが存在しないため、データの著者によって解析されたデータセットを選択した。解析結果を比較することで、異なるツールの性能についていくつかの結論を導き出すことができる。これらの実験では、評価指標としてF1スコア、precision、recall、Jensen-Shannon divergenceを使用した。StrainScanのベンチマークは、Krakenuniq（V0.5.8）[30]、StrainSeeker（V1.5）[31]、Pathoscope2（V2.0.6）[39]、Sigma（V1.0.1）[38]、StrainGE（V1.1.5）[28]、StrainEst（V1.2.4）[29]などの一般的な参照ベースのひずみレベル解析ツールと比較しました。

StrainGEとStrainEstの2つの分解能レベルでの評価これらのテストツールのうち、StrainGEとStrainEstは菌株をクラスターにグループ化し、各クラスターに対して代表的な菌株のみを保持する[28, 29]。したがって、菌株レベルとクラスターレベルの2つの分解能レベルで性能を評価した。菌株レベルの評価では、同定された代表菌株が現在の菌株と同一である場合にのみ、出力を真陽性（TP）としてカウントする。クラスター・レベルの評価では、返された代表株が対象株と同じクラスターにある場合に、出力をTPとしてカウントする。これに対応して，FP の定義もクラスター・レベルではより緩やかになる．他のすべてのツールについては、関連する統計量を計算するために菌株レベルの分解能を使用した。以下に実験結果を示す。

参照データベースの構築
本研究でテストしたすべての生物種について、可能な限り包括的に参照株ゲノムデータベースを作成した。そのため、NCBI RefSeqデータベースから全ゲノムおよびドラフトゲノムをダウンロードした。しかし、25,349の大腸菌ゲノムがあり、1TB以上のメモリを必要とした。しかし、大腸菌のゲノムは25,349個もあり、1TB以上のメモリを必要とする。ハードウェアリソースの制約から、RefSeqの大腸菌完全ゲノムのみを使用した。大腸菌と同様に、結核菌のドラフトゲノムおよび完全ゲノムをすべて使用することは、ハードウェアリソースの制約上不可能である。さらに、入手可能な結核菌ゲノムの中には、10位以下の違いしかないものもある[37]。このようなほぼ同一の菌株は、前処理の段階でクラスタ化される。そこで、Dashing [40]を用いて全結核菌株のペアワイズJaccard類似度を計算し、k-mers Jaccard類似度閾値99％を用いて完全連鎖クラスタリングを行った。そして、他の全ゲノムとの平均類似度が最も高い株のみをそのクラスターに保持した。その結果、6,752ゲノムのうち792ゲノムが結核菌ゲノムとして保存された。

表2 6つの検査対象菌のリファレンスゲノムの要約統計量。「平均Jaccard類似度」は、Dashing[40]を用いて全菌株のk-mers Jaccard類似度の平均を計算することにより得られる。
フルサイズの表
最終的な菌株の数とその他の特性を表2に記録した。すべてのテストツールの入力として使用されたゲノムの数は、"# of input genomes "の列に示されている。前述したように、StrainEstとStrainGEは入力された菌株ゲノムをクラスタリングし、各クラスタから選択された代表的な1株のみを最終的なデータベースに保存する。その結果、残る菌株の数はかなり少なくなります（表2）。StrainEstとStrainGEのクラスターを詳しく見てみると、代表株と同じクラスター内の他の株との間で、遺伝子内容やSNVに大きな違いがあることが観察される（補足図S3、S4、S5）。種によっては、実際の株と代表株との間に5000以上のSNVが存在する。同じクラスター内でも、最も長い株は最も短い株よりも1000個以上の遺伝子を持つことがある。ある最近の研究[56]が示すように、特定の株で見られる「シングルトン」（ユニークな遺伝子）は、株の特性を理解する上で非常に重要である。従って、同じクラスター内の菌株間で、このような遺伝子含有量の大きな差が、異なる性質や機能につながる可能性がある。顕著な例として、病原性株である大腸菌CFT073とプロバイオティクス株である大腸菌Nissle 1917という非常によく似た2つの株が、それぞれStrainEstとStrainGEによって同じクラスターにグループ化されていることが挙げられる。

また、StrainScanでは、クラスタ内の菌株同定を行う前に、菌株をクラスタにグループ分けしている。我々の実験結果から、CSTを用いたクラスタ探索は、テストしたすべての細菌で100％の精度を達成できることが示された。ほとんどの細菌において、StrainScanはStrainEstやStrainGEよりもクラスターの粒度が細かく（補足図S6）、クラスター・レベルでの分解能が高いことを示している。

図6
図6
A 異なるシーケンス深度下での「single-strain」シミュレーションデータセットにおける7ツールのF1スコア。B テストした7ツールの実行時間の比較。SigmaとPathoscope2は、対応する細菌のデータベースを構築したり、シミュレートされたリードから菌株を同定するには計算コストがかかりすぎるため、いくつかのデータセットでは値がない。

フルサイズ画像
シミュレーションリードからの参照株の検出
この実験の目的は、StrainScanやその他のツールを使ってサンプル中の現在の菌株を同定することをテストすることです。課題は2つあります。1つ目の課題は、類似性の高い他の菌株から真の菌株を識別することです。2つ目の課題は、深さの浅い株を識別することです。そこで、シーケンス深度の異なる複数のデータセットを作成した。

それぞれの細菌について、ランダムに参照株を選び、そのシミュレートされたショートリードをすべてのツールの入力として使用した。データに関連するバイアスを避けるため、毎回ランダムに菌株を選び、実験を60回繰り返した。P. copriについては、ゲノム数が少ないため、50回しか実験を繰り返さなかった。選択した各株について、異なるシーケンス深度（10倍、5倍、3倍、1倍）のリードをシミュレートした。その結果、合計1400のデータセットが得られた。各データセットについて、ART [57]を用い、次のパラメータでイルミナリードのシミュレーションを行った： -p -l 250 -f depth -m 600 -s 150、ここでdepthは指定したシーケンス深度である。StrainScanおよび他の6つのプログラムを用いて、これらのシミュレートリードから株を同定した。SigmaとPathoscope2は計算コストが高いため、大腸菌、表皮菌、結核菌についてはデータベースを構築できなかった。

各プログラムのF1スコアを図6Aに示す。また、TP、FN、FP、Recall、PrecisionをSupplementary Table S4に示す。これらの選択された菌株の中には、少なくとも1つの他の参照菌株ゲノムと99.5％以上のk-merベースのJaccard類似性を持つものがある。その結果、いくつかのツールはF1スコアが低い。StrainScanは、深さが1倍以上の場合、すべてのデータセットでほぼ完璧なF1スコアを達成した。しかし、類似性の高い菌株を含む種（A. muciniphilaとM. tuberculosis）では、StrainScanのF1スコアは深度が1Xのときに少し低下する。深度が浅く菌株の類似度が非常に高い場合、CSTアルゴリズムはユニークカバレッジが低いためにいくつかの菌株を同定できず、その結果リコールが低下する。それでもなお、StrainScanはこれらの種で最良のF1スコアを示している。現在のところ、StrainScanが許容する最小深度は1Xであり、深度が1Xより低いと性能は急速に低下する。

Pathoscope2とSigmaは比較的高いリコールを達成しているが、その出力には多くのFPがあり、そのため精度は他のツールよりもかなり低くなっている。Krakenuniqは、多くの菌株がユニークなk-merを持つデータセットで高いF1スコアを達成した。しかし、いくつかの細菌が非常に類似している菌株では、Krakenuniqのリコールと精度が低くなる。StrainSeekerとStrainEstもFPの数が多く、精度が低い。さらに、StrainEstは深さが5倍以下の菌株を同定できない。StrainGEは、多くのデータセットにおいて、クラスタレベルの分解能でStrainScanと同等の性能を示す。しかし、シーケンス深度が浅くなると、より多くのFPを返します。例えば、C. acnesにおけるStrainGEのクラスターレベルの精度は、株深度が10倍から1倍に減少すると、0.95から0.76に低下します。対照的に、StrainScanでは深さが浅くなるにつれてFPが生成されなくなる（補足表S4）。クラスター・レベルでも，結核菌ではk-mersのJaccard類似度が高いため，StrainGEの性能は理想的とはいえない．テストしたツールのうち、StrainSeekerは同じスコアの複数の菌株を返す傾向がある。これは、StrainEstおよびStrainGEが代表的な菌株を返すのと同様であり、菌株グループ間のより細かい区別は提供されない。その結果、これらのツールは解像度が低いという問題がある。例えば、StrainGEは大腸菌で200程度のクラスターについて代表株を返す（補足図S7）。これらのクラスターに含まれる菌株の遺伝的差異を以前分析した結果（補足図S3）から、解像度は理想的とは言えない。

図6Bに異なるツールの実行時間を示す。StrainScanはM. tuberculosisを除くすべての検査対象菌で効率的であった。結核菌の菌株間でk-merベースのJaccard類似性が高いため、StrainScanはほとんどの菌株を、相当数のk-merを持つ1つの大きなクラスターに割り当てた（補足図S6）。それでもなお、StrainScanは結核菌株の同定において最高の再現性と精度を示した。すべての菌株識別実験は、2.4Ghzの14コアIntel Xeon E5-2680v4CPUと128GBメモリを搭載したHPCCのCentOS 6.8ノードでテストした。すべてのツールに8スレッドを使用した。まとめると、StrainScanは、真の菌株が高い配列類似性を持つピアを持つ場合でも、解像度を犠牲にすることなく高い精度を達成することができる。

シミュレーションデータからの共存株の検出
ヒトに関連する微生物叢は、多くの場合、同じ種の近縁の菌株が複雑に混在していることが示されている[15]。Krakenuniq、StrainSeeker、StrainGE、StrainEst、およびStrainScanの同一種の複数菌株の同定性能を定量的に比較するために、6つの細菌から無作為に選択した2、3、および5菌株を含むシミュレートデータセットを作成した。SigmaとPathoscope2はこれらのデータセットの処理に時間がかかりすぎたため、今回の実験には含めなかった。

菌株間の類似性がツールの性能にどのような影響を与えるかを調べるため、複数の菌株の選択に2つの戦略を用いた。StrainScanのクラスタリングステップにおいて、k-merベースのJaccard類似度が95％以上（近似ANI99.89％に相当）の菌株は同じクラスタにグループ化される。そのため、同じクラスターに属する共存菌株を識別・区別することは、異なるクラスターに属する菌株を識別・区別することよりも困難である。異なる難易度を考慮するため、第一の戦略では異なるクラスターからランダムに菌株を選択し、第二の戦略では同じクラスターから異なる菌株を選択した。それぞれの戦略について、ランダムに2、3、5株（3グループ）を選択し、異なるカバレッジプロファイルを使用してショートリードのシミュレーションを行った： 2株は100倍と10倍、3株は100倍と50倍と10倍、5株は100倍と70倍と50倍と20倍と10倍である。その他のリードシミュレーションのパラメータは「単一株」の実験と同じである。次に、別の菌株群を選んで実験を10回繰り返した。最終的に、各細菌種について、1つ目と2つ目の戦略を用いて、異なる菌株数を含む30セットのデータを作成し、合計60セットのデータを作成した。つまり、6つの細菌種について合計360（60×6）のシミュレーション・データセットが存在することになる。

多系統の場合のクラスターレベルの性能評価 StrainGEまたはStrainEstで定義された同じクラスターからn個の菌株が含まれるサンプルの場合、これらのツールの設計に基づき、1つの代表菌株のみが返されます。この1つの代表ストレインを使用すると、多ストレイン実験のリコールが非常に小さくなります。これを避けるため、返された代表株はn回カウントされ、通常、これらのツールのリコールは1.0になります。我々のクラスタはStrainGEやStrainEstで定義されたクラスタよりも粒度が大きいため、同じクラスタの株からシミュレートされたサンプルはすべてこのケースに属します。

図7
図7
A 「複数株」シミュレーションデータセットに対する5つのツールのF1スコア。タイトルの "cluster "は、CSTアルゴリズムによって生成されたクラスターを指す。細菌種ごとに異なる類似性を持つ2株、3株、5株を含む60セットのシミュレートリードがある。StrainSeekerは結核菌株を同定できないため、関連スコアは0である。 B テストした5つのツールの実行時間比較

フルサイズ画像
ベンチマーク結果異なるツールのF1スコア比較を図7Aに示す。TP、FN、FP、recall、precisionは補足表S5、S6に示す。図7Aに示すように、StrainScanはテストしたすべてのデータセットでほぼ完璧なF1スコアを達成した。StrainGEは、P. copriとE. coliの実験でStrainScanと同じクラスタレベルのF1スコアを示した。StrainGEのクラスター・レベルのF1スコアは、一般に入力菌株の類似度が高くなるにつれて低下し、これは図7Aの左と右のパネルを比較することで観察できる。例えば、A. muciniphilaの5株データセットでのクラスター・レベルの再現性は、0.94から0.42に低下した（補足表S6）。単系統の実験と同様に、結核菌に対するStrainGEのクラスタレベルの性能はまだ低い。StrainEstの出力にはFPが多く、そのため精度とF1スコアが低くなっています。しかし、クラスターレベルの性能を評価する方法のため、StrainEstのクラスターレベルのリコールは、同じクラスターの菌株の方が異なるクラスターの菌株よりも高くなり、F1スコアが向上する（図7A右図）。残りのツールは全般的にF1スコアが低い。その中で、Krakenuniqは、同じクラスターに属する菌株の同定よりも、異なるクラスターに属する菌株の同定の方が優れており、これはその手法に沿ったものであった。StrainSeekerは、すべてのテストデータセットでFPが多く、ほとんどのテスト細菌でリコールも低いことから、複数の菌株の同定には不向きであることが示された。さらに、同じスコアで返された菌株の数（StrainSeeker）および返されたクラスター内の菌株の数（StrainEstおよびStrainGE）を分析することにより、菌株間の類似性が高いと、StrainSeeker、StrainGEおよびStrainEstの分解能がさらに低下することもわかった（補足図S8）。全体として、他のテストツールと比較して、StrainScanは高分解能を維持したまま、すべてのデータセットについて系統レベルでのF1スコアで20％以上の改善を達成した。

StrainScanは、結核菌を除き、StrainGEやStrainEstよりも高速である（図7B）。前節で述べたように、結核菌は菌株間の類似性が高いため、StrainScanは同じ大きなクラスター内の菌株を区別するために計算効率を犠牲にしており、その結果、結核菌の処理に時間がかかっている。

図8
図8
グランドトゥルースと予測相対存在量の間の5つのツールのJensen-Shannon Divergence（JSD）。StrainGEとStrainEstは、同じクラスターからの菌株について常に1つの代表菌株を返すため、ここに示す結果は異なるクラスターからの菌株のみを含む。同じクラスターおよび異なるクラスターの株を含む完全な結果は、補足図S9にある。

フルサイズ画像
相対存在量の計算合成データセットにおける予測された菌株プロファイルの精度を測定するために、実際の頻度と推定された頻度との間のJensen-Shannon divergence（JSD）を計算した。予測された相対アバンダンスと真の相対アバンダンスの次元が異なる可能性がある場合は、次元の低い方にゼロを加えてJSDを計算します。StrainGEおよびStrainEstの場合、存在量はクラスター単位で計算される。例えば、1つのサンプルに2つの株があり、その2つの株の代表株が同じであれば、代表株とその存在量を2回使用してJSDを計算します。これにより、片方の菌株の存在量をゼロとするよりもJSDが小さくなる。その結果を図8に示す。すべてのケースで、StrainScanによって再構築された菌株分布は高精度で、Jensen-Shannon divergence (JSD)の中央値は<0.01であった。StrainGEとStrainEstは、菌株数が増えるにつれて組成を定量化する性能が低下した。大腸菌、M. tuberculosis、S. epidermidisのように類似性の高い菌株や多数の菌株を含む種では、StrainScanがStrainGEやStrainEstよりも明らかに優れている。StrainSeekerとKrakenuniqでは、ほとんどのケースでJSDが0.1より大きく、これらのツールでは菌株の構成を正確に定量化できないことを示している。菌株数と類似度が増加するにつれて、ほとんどのテストツールのJSD値の中央値は増加するが、StrainScanのJSDはすべてのケースであまり変動しない。これらの結果は、StrainScanが、類似性の高い菌株を含むサンプルであっても、他のツールよりも複雑なサンプルの組成をよりよく定量できることを示している。

複数菌株の低深度実験複数菌株を低深度で同定する際の各ツールの能力を評価するため、異なるカバレッジプロファイルを使用して、事前に選択した2菌株から追加のショートリードをシミュレートした。次に、これらのデータセットを用いてすべてのツールのベンチマークを行った。その結果、StrainScanは、類似性の高い多くの菌株から低深度の菌株を同定する際に明らかな優位性を示した。例えば、カバレッジ10倍と1倍のC. acnes株を同定した場合、StrainScanは0.98のF1スコアを達成したが、2位のKrakenuniqと3位のStrainGEはそれぞれ0.88と0.81のF1スコアしか達成できなかった（補足表S7参照）。全体として、StrainScanは低深度の複数菌株の同定において競争力を示した（補足表S7およびS8）。この実験に関するその他の詳細は、補足セクション2.1に記載されている。

図9
図9
90の模擬大腸菌データセットと30の実大腸菌データセットにおける、同定された菌株とグランドトゥルース間のマッシュ距離。模擬データセットは、StrainScanのデータベース内の全菌株の中で、実際の菌株と最もよく一致する菌株との間の「1-Mash距離」に従って4つのグループに分けられている。

フルサイズ画像
実際の菌株が参照データベースに含まれていない場合は？
対象菌株がデータベースに含まれていない場合、菌株同定ツールはデータベース内のベストマッチを返すことが期待されます。ベストマッチを同定するための様々なツールの性能を評価するために、NCBIから2022年にリリースされた90の完全な大腸菌ゲノムをダウンロードした（補足表S9）。すべてのツールの大腸菌リファレンスデータベースは2021年までの完全ゲノムを用いて構築されているため、これらのゲノムはすべて我々の構築したデータベースには含まれていない。次に、90の大腸菌ゲノムから10倍カバレッジのリードをシミュレートし、すべてのツールの入力として使用した。各データセットについて、図9の「真実」は、大腸菌全1,433株（表2）における実際の株とそのベストマッチとのMash距離である。次に、図9の最初の4つのパネルに、各データセットについて、返された菌株と実際の菌株の間のMash距離を記録した。これは単一菌株の実験であるため、Mash距離は返送された菌株（代表菌株を含む）と実際の菌株を用いて計算されていることに注意すべきである。これら90ひずみのうち、54ひずみのクラスターがStrainScanの参照データベースから欠落している。それにもかかわらず、StrainScanはこれら54株すべてについて、データベース中の最も適合するクラスターを正しく同定することができる（補足表S9）。さらに、StrainScanで同定された株は、他のツールよりも実際の株とのマッシュ距離が小さい。StrainGEとStrainEstによって同定された最も近いマッチも、"1-MashD "が0.98から0.99の範囲であれば、"truth "とのMash距離はほぼ同じであるが、"1-MashD "が大きくなると性能が低下する。例えば、実際の菌株とベストマッチの間の「1-MashD」が（0.999, 1）の範囲内にある場合、StrainGEとStrainEstは比較的大きなクラスターから代表的な菌株を返す傾向があります（補足表S9）。その結果、これらのツールは、実際の菌株よりもMash距離が大きい代表菌株を返した。対照的に、StrainScanはクラスター内の菌株同定を達成し、より正確なベストマッチを返す。

よりロバストなテストのために、NCBIにあるドラフトゲノムを持つ30個の実際の大腸菌全ゲノムシーケンスデータをダウンロードした。データに関連するバイアスを避けるため、30個の実際のシーケンスデータセットは3つの異なるプロジェクト（PRJNA509690、PRJNA479542、PRJEB21464）からランダムに選択され、それらのカバレッジは20倍から94倍の範囲である[46, 47]。NCBIによると、これらのドラフトゲノムのアセンブリレベルはすべて「Complete」ではなく「Scaffold」である。したがって、これらのゲノムはすべて、完全ゲノムを用いて構築された大腸菌参照データベースにも登録されていない。さらに、これら30種のドラフトゲノムのRefSeqアクセッションを、構築したデータベースのゲノムと比較したが、一致するものは見つからなかった。次に、これらの実データセットにすべてのツールを適用し、同定された株とドラフトゲノムを比較した（図9の最後のパネル）。その結果、StrainScanは他のツールよりも正確なベストマッチを返した。

データベース内の参照配列を欠く複数の菌株がサンプル内に共存している可能性を考慮し、単一菌株実験から得られた90株の大腸菌ゲノムを用いて、追加の2菌株実験を行った。その結果、データベースに参照ゲノムを持たない複数株の同定において、StrainScanは95％のF1スコアを達成したが、2番目に優れたツールであるStrainGE（クラスターレベル）は78％のF1スコアしか達成できなかった（補足図S10参照）。さらに、StrainScanはテストしたすべてのデータセットで偽陽性の同定を行わず、StrainScanによって同定された株は、他のツールよりもグランドトゥルースとのマッシュ距離が小さかった（補足表S10参照）。この実験に関するその他の詳細は補遺2.2節に記載されている。

スパイクされたメタゲノム配列におけるStrainScanの評価
これまでの実験では主にシミュレートした全ゲノムシーケンスデータまたは実際の全ゲノムシーケンスデータを使用したが、ここでは異なる生物種からのリードを含むメタゲノミックデータでもStrainScanが同じ性能を維持するかどうかを評価する。そのために、スパイクしたメタゲノムデータを用いた実験を行った。具体的には、P. copriとE. coliのシミュレートデータセットと実データを混合し、130のスパイクされたメタゲノムデータを作成した。その結果、StrainScanはこれらのスパイクされたメタゲノミックデータセットでも、シミュレートされた全ゲノムシーケンスデータと同じ結果を示し、複雑なサンプルに対するStrainScanの頑健性を実証した（補足表S11）。この実験に関するさらなる詳細は、補足セクション2.3に記載されています。

モックデータにおけるStrainScanの評価
HMPモックデータこの実験では、Human Microbiome Project [49]の2つのサンプルでStrainScanをテストしました。これらのサンプルには、偶数組成（SRR172902）または時差組成（SRR172903）の21の既知の生物が含まれています。21の生物のうち、3つの細菌（E. coli、C. acnes、S. epidermidis）は、菌株レベルの解析が困難なケースであり、我々はこれらの細菌の参照インデックス構造を確立している。そこで、この3つの細菌についてStrainScanを用いて菌株レベルの解析を行った。与えられたデータ記述によると、これらの2つのデータセットでは各細菌は1株しか持っていない。これは単一株の検出であるが、サンプル中の存在量が低いターゲットもあるため、このテストは困難である。比較のために、Krakenuniq、StrainSeeker、StrainGE、およびStrainEstも使用して、これら2つのデータセット内のこれらの細菌の株を同定しました。これら5つのツールの結果を表3に示す。

表3 HMPプロジェクトの2つの模擬コミュニティの解析結果。"#"：出力中に同一のスコアで複数のヒットがある。「NA"：欠損値。太字：予測された優占株の真実とのマッシュ距離が0.05％未満。StrainSeekerでは、真実とのMash距離が最も小さい株を、予測された優性株とした。括弧内の2つの数字は、それぞれ、複数ヒットの数と、すべてのヒットとグランドトゥルース間の平均マッシュ距離を表す
フルサイズの表
テストしたすべての種について、StrainScanは、真性株と非常に類似している（真性株とのMash距離が0.05%未満）1つの優性株の存在を正しく同定した。StrainScanの他に、StrainSeekerとStrainGEも、真正株と非常に類似した株を返した。しかし、StrainSeekerの出力には同一のスコアでヒットした株が複数含まれることが多く、正確な評価が困難である。例えば、StrainSeekerは2つのテストデータセットで119株の大腸菌を返したが、これではユーザーはこれらのサンプルに含まれる実際の菌株を知ることが難しい。表3では、StrainSeekerが予測した優勢株を、真実とのMash距離が最小の株としています。StrainGEもまた、真理値とのMash距離が小さい菌株を返します。残りの2つのツールのうち、StrainEstは低存在株を同定できず、実行に時間がかかった。

複数株を含む大腸菌モックコミュニティ複数株の構成が既知の実サンプルに対する各ツールの性能を評価するため、宿主（すなわちヒト）と4つの異なる大腸菌株からの大量のリードを含むモックコミュニティシーケンスデータセット（SRR13355226）をダウンロードした。すべてのリードはすべてのツールの入力として使用される。StrainScanは、偽陽性の同定がなく、菌株レベルで4つの菌株を正しく同定した唯一のツールであった（補足図S11）。StrainGEは代表的な4株を正しく同定し、偽陽性の同定はなかった。しかし、StrainGEによって同定された代表株と真の株との間には136の異なる遺伝子が存在する。これらの遺伝子は Prokka [58]を用いて予測され、比較解析は Roary [59]によって終了した。StrainEstも4つの代表株を正しく同定したが、多くのFPを返した。残りの2つのツールは、4つの菌株をすべて正しく同定することはできず、中でもStrainSeekerの出力には同一のスコアで複数のヒットが含まれ、Krakenuniqは多くの偽陽性菌株を報告した。

StrainScanは実際のシーケンスデータから病原性株を検出する
病原体検出におけるStrainScanの応用の可能性を示すため、StrainScanを適用して、2つの研究（BioProject Accession: PRJEB1775およびPRJEB2777）における大腸菌と結核菌の病原性株の存在を調べた。最初の研究は、2011年にドイツで発生した大腸菌アウトブレイク[50]に関連するもので、腸内凝集性（EAEC）株が原因であった。これらのデータセットは、便サンプルからメタゲノム配列決定を用いて配列決定されている。原著論文によると、病原株は大腸菌O104:H4である。2つ目の研究は、未治療の結核患者に対する治療を評価したREMoxTB臨床試験の患者における結核の再発頻度を調査したものである[51]。シーケンスデータは細菌分離株から取得され、各データセットには、サンプルに含まれる菌株を表すアセンブルされたゲノムが公開されています。これら2つの研究それぞれから、実験用に6つのサンプルを選択した。比較のため、これらの実シーケンスデータから病原性株を検出する他のツールも適用した。表4に示すように、StrainScanはテストしたすべてのデータセットで正しい菌株を同定できたが、他のツールは一部のデータセットで正しい菌株を検出できなかった。StrainEstはほとんどのデータセットで正しい菌株を同定できたが、一部のデータセットでは正しい菌株の代表的な菌株しか返さなかった。

表4 大腸菌と結核菌の病原性株の同定におけるStrainScan、StrainGE、KrakenUniq、StrainEst、StrainSeekerの性能。緑：グランドトゥルースと一致した結果。赤: グラウンドトゥルースと一致しない結果。"#"：出力に同じスコアを持つ複数のヒットがある。"R"：同定された菌株がグランドトゥルースの代表菌株である。StrainSeekerは結核菌株を同定できないため、関連する結果は"-"であることに注意。StrainSeekerの場合、括弧内の数字は複数ヒットの数を表す。
フルサイズの表
StrainScanは抗生物質耐性S. epidermidis株の動態を正確に検出する
この実験では、公開メタゲノムデータセット（PRJNA490375）中の2つのS. epidermidis株の動態を検出するStrainScanの能力をテストした。オリジナルの研究[52]によると、これらのデータセットはin vitroでの混合株培養から作成されている。著者らはS. epidermidisの2つの皮膚分離株を1:1の比率で2つのグループに分けて培養した。一方のグループは抗生物質エリスロマイシン処理（Ery）で培養し、もう一方は抗生物質処理なしで培養した（no_ATB）。これら2つの分離株のうち、NIHLM023は抗生物質エリスロマイシンに耐性を示さなかったが、NIHLM001はエリスロマイシンに強い耐性を示した。最後に、3つの異なる時点でのメタゲノムシークエンシングにより、2つのグループから6つのデータセットが得られた。これら2つの菌株の相対的な存在量は示されていないが、各時点でのカバー率が示されており、これを用いて2つの菌株の割合の変化を評価することができる。原著論文で報告されたカバレッジ比によると、各時点でno_ATB群ではNIHLM023株が常に優勢であったのに対し、Ery群ではNIHLM001株が優勢であった。次に、これら6つのデータセットにStrainScanと他のツールを適用した結果を図10Aに示す。StrainScanは、各時点で2株の正しい株比率を返し、すべてのサンプルで偽陽性の同定がない唯一のツールである。テストしたツールのうち、StrainGEとStrainEstは一部のサンプルで正しい代表株を返した。そこで、これらの代表株と実際の株との間で異なる遺伝子を調べた。遺伝子はProkka[58]で予測し、比較解析はRoary[59]で行った。図10BおよびCに示すように、これらの菌株の間にはまだ多くの異なる遺伝子が存在する。これらの差異遺伝子の中には、菌株の機能にとって非常に重要なものもある。例えば、NIHLM001はssaA_1という遺伝子を持っており、この遺伝子は薬剤耐性などこの菌株の多くの性質と関連していることが示されている[60,61,62]。しかし、この遺伝子は代表的な菌株であるUMB8493には備わっていない。より正確な菌株レベルの解析には、菌株レベルの分解能が望ましい。

図10
図10
A 6つの実メタゲノムサンプルにおける5つのツールによるS. epidermidis株の推定存在量。No_ATB：抗生物質無投与群、非耐性株NIHLM023が優勢株。Ery：抗生物質エリスロマイシン処理群で、耐性株NIHLM001が優占株。各色は菌株を表す。B-C StrainGEおよびStrainEstで同定された代表株と実際の株との異なる遺伝子数。

フルサイズ画像
StrainScanはメタゲノムデータから低深度病原性C. difficile株を同定する
低深度病原性株を同定するStrainScanの能力をさらに検証するために、低深度病原性Clostridioides difficile株を含む1つのメタゲノミックデータセットにStrainScanと他のツールを適用した[53]。オリジナルの研究によると、このデータセットからは深さの浅い2つの病原性C. difficile株（ \sim1X ）が検出された。これら2つの菌株の間には3つの変異しか報告されていない。別の研究[63]でもC. difficile株間で同じ変異が検出されており、配列決定エラーというよりも、類似性の高い複数の株が存在した結果であることが示唆された。C. difficileは6つの対象菌の1つではないため、まずNCBI RefSeqからダウンロードした102の完全なClostridioides difficileゲノムを用いて、各ツールの参照データベースを構築した。StrainScan、StrainGE、Krakenuniq、Pathoscope2、Sigmaは約1倍の深さで菌株を検出できた（表5）。しかし、Krakenuniq、Pathoscope2、Sigmaは10株以上を同定しており、多数の偽陽性が存在することを示している。StrainGEは1株しか出力しなかったが、同定された株は病原性Clostridioides difficileが持つはずの2つの遺伝子TcdAとTcdBを含んでいなかった。StrainScanは、偽陽性を示すことなく低深度の強毒株を検出した唯一のツールであった。しかし、StrainScanもStrainGEも、深度が浅く、優勢株との類似性が非常に高いため、もう一方の菌株を見逃していた。

表5 低深度病原性Clostridioides difficile株を含む1つのメタゲノムデータセットに対する7ツールの同定結果。「NA」：欠損値。「TcdA」と「TcdB」：病原性C. difficile株の2つの重要な遺伝子。TcdA "と "TcdB "の存在はCarbohydrate-Active enzymes database [64]を用いて検証されている。
原寸表
図11
図11
StrainScanにより、9つの実際のメタゲノミックサンプルにおけるC. acnesの多様性が明らかになった。これらのサンプルは、2人の健常人（HV05とHV06）の皮膚の3つの異なる部位（AI、Lc、Vf）から、異なる時点で採取された。部位コードは原著[15]に記載されている。B 同定された菌株の系統樹。葉は(A)と同じスキーマで色付けされ、距離はMash距離[14]である。ツリーはiTOL [65]によって可視化されている。C "GCF_000194905 "と "GCF_000252385 "のユニークな遺伝子クラスターの存在。「GCF_000194905 "と "GCF_000252385 "は類似性が高く、StrainGEとStrainEstの同じクラスターに存在する。存在はオリジナルの研究[15]によって与えられ、各列は1つのユニークな遺伝子クラスターを指す。

フルサイズ画像
StrainScanにより、ヒト皮膚におけるC. acnesの多様性が明らかになった。
先行研究[11, 15]によると、C. acnesはヒトの皮膚上で最も一般的な細菌の一つであり、通常、複雑な複数菌株群集を有する。StrainEstは、研究[29]のヒト皮膚データセット（SRP002480）の再解析に適用された。しかし、一部のサンプルでは、StrainEstで同定された菌株数が元の研究で報告された菌株数より少ないことが多いことが判明した。オリジナルの研究では、著者らは社内のパイプラインを使用して菌株を決定し、サンプル中の相対的な存在量を予測しており、この社内のパイプラインの精度は、ヒトの皮膚マイクロバイオームデータを対象とした広範なシミュレーションで検証済みであった[11, 15]。そこで、このケースに当てはまる2人の個人から9つのサンプルを選択し、StrainEstとStrainScanを使用してこれらのサンプルを再分析した。そして、それらの出力を元の研究で報告されたものと比較した。StrainGEの競争力を考慮して、StrainGEも比較に加えた。一貫性を保つため、StrainEst、StrainGE、およびStrainScan用の新しいカスタムデータベースを構築するために、元の研究で使用されたすべての株を選択し、以降の解析にはこれらの新しく構築されたデータベースを使用した。その結果を図11に示す。StrainEstおよびStrainGEでは、すべてのテストデータセットで1～2株しか返されなかったが、原著論文ではより多くの株が報告されている。StrainScanはStrainEstやStrainGEよりも、報告された結果と類似した相対存在量パターンを示した。また、StrainScanでは、「HV05_AI」のサンプル中に、元の研究では検出されなかった菌株も検出され、C. acnesの多様性がより高いことが示唆された。例えば、「GCF_000145115」と「GCF_004136215」という2つの類似性の高い菌株が、それぞれStrainScanとStrainGE/StrainEstによって同定された。図11Bのグレーのボックスで示すように、両者はツリー上で非常に近接しており、"GCF_004136215 "はStrainGEおよびStrainEstにおいて "GCF_000145115 "の代表株である。しかし、両者は76の異なる遺伝子を持つ。原著論文[11]の SNP 解析によると，"2011-12-14 "のサンプル "HV05_AI "では，"GCF_000145115 "の 22 個のユニーク SNP が同定されたが，"GCF_004136215 "のユニーク SNP は検出されなかった．したがって、このサンプルには "GCF_004136215 "ではなく、"GCF_000145115 "が存在する可能性が高く、StrainScanの高分解能の優位性が示された。

また、同一サンプル "HV06 "に対して、StrainGEとStrainEstは異なる代表株を同定した（図11A）。この結果は、クラスタベースの手法の同定結果が、クラスタリングや代表株の選択戦略の違いによって影響を受ける可能性があることを示している。一方、StrainScanでは、これらの試料から "GCF_000194905 "と "GCF_000252385 "という2つの類似性の高い菌株が同定されており、報告された結果と一致している。異なる菌株構成の結果は、下流の分析に影響を与える可能性がある。ここでは、"HV06_Vf "の解析結果を一例として取り上げた。GCF_000194905 "と "GCF_000252385 "のユニークな遺伝子クラスターにリードをアライメントし、"HV06_Vf "に含まれるこれら2株の遺伝子含量の経時変化を解析した。図11Cに示すように、2株由来のユニークな遺伝子クラスターの存在は、異なる時点で変化しており、この個体における株レベルの機能変化を反映していた。この結果は、StrainScanが類似性の高い系統を区別することを示しており、これらの系統を一緒に用いることで、より包括的な系統レベルの機能解析を行うことができる。

実際のメタゲノムシーケンスデータにおけるStrainScanのさらなる応用を示す2つの例
実際のメタゲノムシークエンシングデータにおけるStrainScanの幅広い応用例を示すため、横断的研究における大腸菌株の解析[48, 54]と、異なる集団におけるP. copri株の解析[9]にStrainScanを適用した。解析結果から、StrainScanはメタゲノミックサンプルからより高い解像度で正確に菌株を同定でき、より包括的な生物学的洞察につながることが示された。例えば、StrainScanは3カ国のサンプル中の大腸菌株を3つの異なるグループに区別できる唯一のツールである（補足図S12）。同様に、StrainScanで同定されたP. copri株を分析したところ、雑食性（Omnivores）サンプルの株と菜食性（Vegans）サンプルの株は系統関係の点で明確に分離しており、菜食性（Vegetarians）サンプルの株はその中間に位置していることが観察された（補足図S12）。これら2つの実験に関する詳細は、補足セクション2.4および2.5に記載されている。

考察
本研究では、ショートリード用の新しい菌株レベル組成解析ツールであるStrainScanを紹介した。系統同定精度と計算の複雑さのバランスをとるために、ツリーベースのk-mersインデックス構造を新たに設計した。そして、情報量の多いk-mersと弾性ネットモデルを適用して、ひずみを同定し、その存在量を予測することで、StrainScanはひずみレベル解析の分解能と存在量推定の精度を向上させた。

StrainScanは、複雑さの異なるすべてのベンチマークデータセットにおいて、他のテストツールよりも高い精度と解像度を示した。特に、StrainScanは、より高い類似性と様々な配列決定深度の株を含むデータセットにおいて、他の全てのツールよりも優れている。このレベルの高分解能は、Pathoscope2やSigmaのようなアライメントベースのツールでも達成できる。しかし、StrainScanはそれらよりも少なくとも10倍速い。模擬データと実データの実験結果は、StrainScanがより包括的な株レベルの組成分析を提供できることをさらに実証している。

参照データベースにない新規の菌株はしばしば出現する傾向がある。StrainGEやStrainEstのようなツールは、サンプル中に存在する菌株に最も類似した代表的な菌株を返しますが、StrainScanは最もよく一致する菌株を返します。その結果、StrainScanが返す菌株はより高い精度で実際の菌株を表すことができる。さらに、StrainScanは同定されたクラスターに関する情報を出力するため、実際の菌株のゲノム配列が入手できない場合でも、同定されたクラスターを下流の解析対象として検討することができる。ベンチマーク実験（図9および補足図S10）に示されるように、StrainScanは、サンプル中の実際の菌株とのマッシュ距離がより小さい菌株を同定することができ、下流の解析により正確な参照菌株を提供することができる。その上、図11Cに示すように、類似性の高い参照株は、依然として株特異的遺伝子を含んでいる可能性があるため、（代表的な株ではなく）すべての株を使用することで、複数のサンプルの比較において、より包括的な株レベルの機能変化を明らかにすることができる。これらの結果は、StrainScanによって同定された菌株を下流解析に用いることで、代表的な菌株のみを返すツールと比較して、新たな生物学的洞察を生み出す大きな可能性があることを示している。

菌株が高度に類似している場合、シーケンスの深さがStrainScanの分解能に影響する主な要因である。したがって、菌株の類似度が高く、配列深度が浅い場合、CSTを用いたクラスター探索では、カバレッジが低いためにいくつかのクラスターを同定できない。その結果、菌株の深さが1倍を下回ると、StrainScanの性能は急激に低下する。さらに、StrainScanが必要とする最小配列決定深度と、低深度サンプルで株を検出する可能性を調査した（補足セクション2.6）。25,349の完全およびドラフト大腸菌ゲノムのCSTを構築して、StrainScanの限界をテストした。このプログラムには1TB以上のメモリが必要である。我々の経験的テストから、StrainScanは5000ゲノム以下のCST構築に効率的であることが示された（補足表S12参照）。しかし、ゲノム数が多く、種内の多様性が高い細菌（例えば、大腸菌）では、CSTの構築に時間がかかることがあります。そのため、5000以上のゲノムがある場合は、完全なゲノムのみを用いてインデックス構造を構築することを推奨する。クラスタ検索に基づく最初のステップは、検索空間を効率的に縮小することができる。しかし、すべての参照ゲノムが非常に類似しており、ほんのわずかな塩基の違いしかない場合、それらは1つのクラスターにまとめられる傾向がある。このような場合、クラスター検索は依然として大きな検索空間を第2ステップに返し、クラスター索引構造を十分に活用できない。M. tuberculosisはクラスターが大きいため、菌株同定が遅くなる。参照ゲノムの質は、StrainScanを含む参照ベースの手法[43]に影響を与える可能性があることに注意すべきである。StrainScanでは、ユーザーが自分のゲノムを用いてデータベースを構築することができます。したがって、入力ゲノムのデータ汚染やバイアスの問題を軽減するための前処理を行うことができる。

上述したように、入力参照ゲノムのスケーラビリティはStrainScanの限界となり得る。しかし、マイクロバイオーム株レベル解析の拡張性を高めるために、多くの効率的なデータ構造[66]や手法[12]が開発されている。例えば、mSWEEP [67]と呼ばれるツールは、擬似アライメント技術を利用することで、大規模な参照ゲノムから菌株系統を正確に同定することができる。したがって、今後の課題として、擬似アライメントに基づく手法やHyperLogLog [30]のような効率的なデータ構造を拡張してStrainScanのスケーラビリティを向上させ、StrainScanのデータベース構築と同定プロセスの両方を高速化することを計画している。

結論
結論として、我々は正確で効率的、かつ高分解能のショートリード用系統レベル組成解析ツールを開発した。シミュレートデータセットとモックデータセットを用いた実験結果から、StrainScanは高分解能を維持したまま、既存のツールよりも正確な菌株レベルのマイクロバイオーム組成解析を実現できることが示された。また、これらのベンチマーク実験の結果は、複雑なサンプル、低存在株、参照データベースにない株に対するStrainScanの頑健性も証明している。さらに、実際のメタゲノミックデータセットの解析結果は、StrainScanが病原性菌株の同定、菌株レベルの組成・機能解析、異なる研究やデータセットにまたがるメタ解析に役立つことを示している。これらの結果から、StrainScanはこの分野への重要な貢献であり、最先端のツールよりも性能が向上していると考えている。

データと資料の入手
本研究に使用したすべてのデータとコードはオンラインで入手可能である。StrainScanのソースコードはhttps://github.com/liaoherui/StrainScan。

参考文献
Luo C, Walk ST, Gordon DM, Feldgarden M, Tiedje JM, Konstantinidis KT. 環境大腸菌のゲノム配列決定により、モデル細菌種の生態と種分化の理解が拡大。Proc Natl Acad Sci U S A. 2011;108(17):7200-5.

論文

PubMed

パブメッドセントラル

Google Scholar

この研究により、野生プロクロロコッカスには数百のサブ集団が共存していることが明らかになった。Science. 2014;344(6182):416–20.

論文

パブコメ

グーグル

ヒト腸内細菌叢のゲノム変異ランドスケープ。Nature. 2013;493(7430):45–50.

論文

パブコメ

グーグル奨学生

生後数カ月間の母子間細菌伝播の菌株レベルでの解析（Yassour M, Jason E, Hogstrom LJ, Arthur TD, Tripathi S, Siljander H, Selvenius J, Oikarinen S, Hyöty H, Virtanen SM, et al. Cell Host Microbe. 2018;24(1):146-54.

論文

CAS

PubMed

パブメドセントラル

Google Scholar

Rasko DA, Webster DR, Sahl JW, Bashir A, Boisen N, Scheutz F, Paxinos EE, Sebra R, Chin CS, Iliopoulos D, et al. ドイツで発生した溶血性尿毒症症候群の原因となった大腸菌の起源。N Engl J Med. 2011;365(8):709-17.

Nayfach S, Rodriguez-Mueller B, Garud N, Pollard KS. 菌株プロファイリングのための統合メタゲノミクスパイプラインにより、細菌の伝播と生物地理学の新規パターンが明らかになった。Genome Res.

論文

PubMed

パブメドセントラル

Google Scholar

Mills RH, Vázquez-Baeza Y, Zhu Q, Jiang L, Gaffney J, Humphrey G, Smarr L, Knight R, Gonzalez DJ. クローン病患者の4.5年間の研究におけるメタプロテオームのメタゲノム予測の評価。2019;4(1):e00337-18.

Tett A, Huang KD, Asnicar F, Fehlner-Peach H, Pasolli E, Karcher N, Armanini F, Manghi P, Bonham K, Zolfo M, et al. Prevotella copri complexは、欧米化した集団では十分に代表されない4つの異なるクレードからなる。Cell Host Microbe. 2019;26(5):666-79.

論文

CAS

PubMed

パブメドセントラル

Google Scholar

ヒト腸内プレボテラ属コプリ株の遺伝的および機能的形質は、異なる習慣食と関連している。Cell Host Microbe. 2019;25(3):444-53.

論文

PubMed

グーグル奨学生

Zhai R, Xue X, Zhang L, Yang X, Zhao L, Zhang C. Strain-specific anti-inflammatory properties of two Akkermansia muciniphila strains on chronic colitis in mice. Front Cell Infect Microbiol. 2019;9:239.

論文

CAS

PubMed

パブメドセントラル

グーグル奨学生

ヒトの皮膚メタゲノムにおける生物地理学的機能と個体差の形成。ヒト皮膚メタゲノムにおける生物地理学と個体性の機能形成。Nature. 2014;514(7520):59–64.

論文

PubMed

パブメドセントラル

グーグル

Schaeffer、Pimentel H、Bray N、他。メタゲノムリード割り当てのためのシュードアライメント。Bioinformatics. 2017;33(14):2082-8.

論文

PubMed

パブメドセントラル

Google Scholar

Sharon I, Morowitz MJ, Thomas BC, et al. 時系列コミュニティゲノム解析により、乳児の腸内コロニー形成期における細菌種、菌株、ファージの急速なシフトが明らかになった。Genome Res. 2013;23(1):111-20.

論文

PubMed

パブメドセントラル

グーグル奨学生

ゲノムとメタゲノム間の距離の推定にMinHashを使用。Mash：MinHashを用いたゲノムおよびメタゲノム距離の高速推定。Genome Biol.

論文

PubMed

パブメドセントラル

Google Scholar

ヒト皮膚マイクロバイオームの時間的安定性。ヒト皮膚マイクロバイオームの時間的安定性。Cell. 2016;165(4):854-66.

論文

CAS

PubMed

パブメドセントラル

Google Scholar

Costea PI, Coelho LP, Sunagawa S, et al. Subspecies in the global human gut microbiome. Mol Syst Biol.

論文

PubMed

パブメドセントラル

Google Scholar

Luo C, Knight R, Siljander H, Knip M, Xavier RJ, Gevers D. ConStrainsはメタゲノムデータセット中の微生物株を同定する。Nat Biotechnol. 2015;33(10):1045-52.

論文

CAS

PubMed

パブメドセントラル

Google Scholar

Truong DT, Tett A, Pasolli E, Huttenhower C, Segata N. Microbial strain-level population structure and genetic diversity from metagenomes. Genome Res. 2017;27(4):626-38.

論文

CAS

PubMed

パブメドセントラル

Google Scholar

MetaMLST: multi-locus strain-level bacterial typing from metagenomic samples. Nucleic Acids Res.

論文

Google Scholar

metaSNV: A tool for metagenomic strain level analysis. PLoS ONE. 2017;12(7):0182392.

論文

Google Scholar

Scholz M, Ward DV, Pasolli E, Tolio T, Zolfo M, Asnicar F, Truong DT, Tett A, Morrow AL, Segata N. ショットガンメタゲノミクスによる菌株レベルの微生物疫学と集団ゲノミクス。Nat Methods. 2016;13(5):435-8.

論文

CAS

PubMed

Google Scholar

微生物群集とヒトマイクロバイオームの菌株レベルの疫学。Genome Med. 2020;12(1):71.

論文

PubMed

パブメドセントラル

Google Scholar

コロニーおよびメタゲノム配列データにおける菌株レベルの微生物検出のための計算手法。Front Microbiol. 2020;11:1925.

論文

パブメドセントラル

Google Scholar

大腸菌の系統学的に近縁なプロバイオティック株と病原株との間に見られる細胞外代謝の大きな違い。Front Microbiol. 2019;10:252.

論文

グーグル・スカラー

Piel D, Bruto M, Labreuche Y, et al. 自然集団におけるファージ-宿主共進化。Nat Microbiol. 2022;7(7):1075-86.

論文

CAS

パブコメ

Google Scholar

微生物種内の多様性：マイクロバイオームにおける菌株の解釈。Nat Rev Microbiol. 2020;18(9):491-506.

論文

PubMed

パブメドセントラル

Google Scholar

ベニーS、ロドリゴDAT、マハルジャンRP、トーマスF. 研究室間での細菌株の移動がもたらす不確実な結果-rpoSの不安定性を例として-. BMC Microbiol. 2011;11:248.

論文

Google Scholar

StrainGE: 複雑な微生物群集における低存在株を追跡し、特徴付けるためのツールキット。Genome Biol.

論文

PubMed

パブメドセントラル

Google Scholar

メタゲノムシークエンシングから得られた細菌種の株プロファイリングと疫学。Nat Commun. 2017;8(1):2260.

論文

PubMed

パブメドセントラル

Google Scholar

Breitwieser FP, Baker DN, Salzberg SL. KrakenUniq：ユニークなk-merカウントを用いた確信のある高速メタゲノム分類。Genome Biol.

論文

CAS

PubMed

パブメドセントラル

Google Scholar

StrainSeeker: ユーザー提供のガイドツリーを用いてシーケンスリードから細菌株を高速同定。PeerJ. 2017;5:3353.

論文

Google Scholar

Quince C, Delmont TO, Raguideau S, Alneberg J, Darling AE, Collins G, Eren AM. DESMAN：メタゲノムから菌株をde novo抽出するための新しいツール。Genome Biol.

論文

PubMed

PubMed Central

Google Scholar

Pulido-Tamayo S, Sánchez-Rodríguez A, Swings T, Van den Bergh B, Dubey A, Steenackers H, Michiels J, Fostier J, Marchal K. 細菌集団のディープシーケンスデータからの頻度に基づくハプロタイプ再構築。Nucleic Acids Res.

論文

Google Scholar

Koslicki D, Falush D. MetaPalette: a k-mer painting approach for metagenomic taxographic profiling and quantification of novel strain variation. 2016;1(3):e00020-16.

Smillie CS, Sauk J, Gevers D, Friedman J, Sung J, Youngster I, Hohmann EL, Staley C, Khoruts A, Sadowsky MJ, et al. Strain tracking reveals the determinants of bacterial engraftment in the human gut following fecal microbiota transplantation. Cell Host Microbe. 2018;23(2):229-40.

論文

CAS

PubMed

パブメドセントラル

Google Scholar

estMOI：寄生虫のディープシーケンスデータを用いた感染の多重性の推定。Bioinformatics. 2014;30(9):1292-4.

論文

PubMed

パブメドセントラル

Google Scholar

QuantTB - 全ゲノムシークエンシングデータ内で結核菌混合感染を分類する方法。BMC Genomics. 2020;21(1):80.

論文

PubMed

パブメドセントラル

Google Scholar

Sigma: バイオサーベイランスのためのメタゲノム解析からのゲノムの株レベル推論。Bioinformatics. 2015;31(2):170-7.

論文

PubMed

Google Scholar

Hong C, Manimaran S, Shen Y, Perez-Rogers JF, Byrd AL, Castro-Nallar E, Crandall KA, Johnson WE. PathoScope 2.0：環境または臨床シーケンスサンプルにおける菌株同定のための完全な計算フレームワーク。Microbiome. 2014;2:33.

Baker DN, Langmead B. Dashing: Fast and accurate genomic distances with HyperLogLog. Genome Biol.

論文

PubMed

パブメドセントラル

Google Scholar

Marcais G, Kingsford C. A fast, lock-free approach for efficient parallel counting of occurrences of k-mers. Bioinformatics. 2011;27(6):764-70.

論文

パブコメ

パブメドセントラル

Google Scholar

ホールGA, スピードTP, ウッドラフCJ. ロングリードとmapqスコアを用いた菌株レベルのサンプル特性解析。bioRxiv.

VirStrain: A strain identification tool for RNA viruses. Genome Biol.

論文

パブコメ

パブメッドセントラル

グーグル奨学生

SibeliaZを用いたスケーラブルな複数全ゲノムアラインメントと局所的に共線的なブロック構築。Nat Commun. 2020;11(1):6327.

論文

PubMed

パブメドセントラル

グーグル奨学生

Zingali T, Reid CJ, Chapman TA, Gaio D, Djordjevic SP. 母豚とその子豚のクラス1インテグロン保有豚糞便常在性大腸菌の全ゲノム配列解析。Microorganisms. 2020;8(6):843.

Cummins ML, Reid CJ, Chowdhury PR, Bushell RN, Djordjevic SP. クラス1インテグラーゼ遺伝子を保有するオーストラリアの鳥類病原性大腸菌の全ゲノム配列解析。Microb Genom 2019;5(2):e000250.

Reid CJ, Wyrsch ER, Chowdhury PR, Zingali T, Djordjevic SP. 豚の常在性大腸菌：IS26に関連するクラス1インテグロンのリザーバー。Microb Genom. 2017;3(12):e000143.

Qin J, Li Y, Cai Z, Li S, Zhu J, Zhang F, Liang S, Zhang W, Guan Y, Shen D, et al. A metagenome-wide association study of gut microbiota in type 2 diabetes. Nature. 2012;490(7418):55–60.

論文

パブコメ

グーグル奨学生

Huttenhower C, Gevers D, Knight R, Abubucker S, Badger JH, Chinwalla AT, Creasy HH, Earl AM, FitzGerald MG, Fulton RS, et al. 健康なヒトのマイクロバイオームの構造、機能、多様性。Nature. 2012;486(7402):207–14.

論文

CAS

Google Scholar

Loman NJ, Constantinidou C, Christner M, Rohde H, Chan JZ, Quick J, Weir JC, Quince C, Smith GP, Betley JR, et al. 志賀毒素原性大腸菌O104:H4の集団発生を調査するための培養に依存しない配列ベースのメタゲノミクスアプローチ。JAMA. 2013;309(14):1502–10.

論文

パブコメ

Google Scholar

Bryant JM、Harris SR、Parkhill J、Dawson R、Diacon AH、van Helden P、Pym A、Mahayiddin AA、Chuchottaworn C、Sanne IM、他。結核菌の再発または再感染を確定するための全ゲノム配列決定：レトロスペクティブ観察研究。Lancet Respir Med. 2013;1(10):786-92.

論文

CAS

PubMed

パブメドセントラル

Google Scholar

メタゲノム解析によるin situでの菌株の増殖速度の推定。サイエンス・アドバンス2020;6(17):2299.

論文

Google Scholar

Džunková M, Moya A, Chen X, Kelly C, D'Auria G. FACSと超低入力ゲノムシーケンスによる混合株感染の検出。Gut Microbes. 2020;11(3):305-9.

Vatanen T, Kostic AD, d'Hennezel E, Siljander H, Franzosa EA, Yassour M, Kolde R, Vlamakis H, Arthur TD, Hämäläinen A, et al. Microbiome LPS免疫原性のばらつきはヒトの自己免疫に寄与する。Cell. 2016;165(6):1551.

論文

PubMed

Google Scholar

Fuglede B, Topsøe F. Jensen-shannon divergence and hilbert space embedding. 情報理論国際シンポジウム. ISIT 2004. Proceedings. 2004;2004:31.

Google Scholar

Tierney BT, Yang Z, Luber JM, Beaudin M, Wibowo MC, Baek C, Patel CJ, Kostic AD. 腸内および口腔内のヒトマイクロバイオームにおける遺伝的内容のランドスケープ。Cell Host Microbe. 2019;26(2):283-95.

論文

CAS

PubMed

パブメドセントラル

Google Scholar

Huang W, Li L, Myers JR, Marth GT. ART：次世代シーケンサーリードシミュレーター。Bioinformatics. 2012;28(4):593-4.

論文

パブコメ

Google Scholar

Prokka: 迅速な原核生物ゲノムアノテーション。Bioinformatics. 2014;30(14):2068-9.

論文

パブコメ

Google Scholar

原核生物パンゲノム解析のための迅速な大規模解析ツールRoary。Bioinformatics. 2015;31(22):3691-3.

論文

PubMed

パブメドセントラル

Google Scholar

Lang S, Livesley MA. 表皮ブドウ球菌の新規抗原の同定。FEMS Immunol Med Microbiol. 2000;29(3):213-20.

論文

パブコメ

Google Scholar

バイオフィルムおよび浮遊性条件下で培養した黄色ブドウ球菌の異なる遺伝子発現プロファイリング。Appl Environ Microbiol. 2005;71(5):2663-76.

論文

PubMed

パブメドセントラル

Google Scholar

ウディンJ、ダワンJ、チョンG、ユーT、アンJ. 抗生物質耐性の拡散における細菌膜小胞の役割と治療薬デリバリーの有望な担体としての役割。Microorganisms. 2020;8(5):670.

クロストリジウム・ディフィシル（Clostridium difficile）のTcdAとTcdBの配列変異： TcdAを切断したST37は中国における潜在的流行株である。J Clin Microbiol. 2014;52(9):3264-70.

論文

CAS

PubMed

パブメドセントラル

Google Scholar

Vincent L, Hemalatha GR, Elodie D, Coutinho PM, Bernard H. 2013年の糖質活性酵素データベース（CAZy）。Nucleic Acids Res. 2014;42(Database issue):490-5.

Letunic I, Bork P. Interactive Tree Of Life (iTOL) v4: 最近の更新と新たな展開。Nucleic Acids Res. 2019;47(W1):256-9.

論文

グーグル・スカラー

Wandelt S, Starlinger J, Bux M, Leser U. Rcsi： 1000(s)単位のゲノムにおけるスケーラブルな類似性検索。Proc VLDB Endow. 2013;6(13):1534-45.

論文

Google Scholar

Wellcome Open Res.

論文

Google Scholar

参考文献のダウンロード

謝辞
該当なし。

資金提供
本研究は、香港城市大学（Project 9678241）、香港革新技術委員会（InnoHK Project CIMDA）、香港研究助成委員会（RGC）一般研究基金（General Research Fund）11206819の支援を受けている。

著者情報
著者ノート
Herui LiaoとYongxin Jiは同等に本研究に貢献した。

著者および所属
香港城市大学電気工学科、九龍、中国

Herui Liao, Yongxin Ji & Yanni Sun

貢献
YSがアイデアを考案し、研究を監督した。HLとYJはStrainScanパッケージの設計と実装を行った。HLは実験を行い、Methods and Resultsの第一稿を執筆した。YJはMethodsの一部を執筆した。YS、HL、YJは原稿を修正した。最終原稿は著者全員が読み、承認した。

責任著者
Yanni Sunまで。

倫理申告
倫理承認および参加同意
該当なし。

発表の同意
すべての著者が掲載に同意している。

競合利益
著者らは、競合する利益はないと宣言している。

追加情報
出版社ノート
シュプリンガー・ネイチャーは、出版された地図の管轄権の主張および所属機関に関して中立を保つ。

補足情報
追加ファイル1.
補足セクション、図、表。補足情報は追加PDFファイルに含まれています。

権利と許可
オープンアクセス本論文は、クリエイティブ・コモンズ表示4.0国際ライセンスの下でライセンスされている。このライセンスは、原著者および出典に適切なクレジットを付与し、クリエイティブ・コモンズ・ライセンスへのリンクを提供し、変更が加えられた場合にその旨を示す限り、いかなる媒体または形式においても、使用、共有、翻案、配布、複製を許可するものである。この記事に掲載されている画像やその他の第三者の素材は、その素材へのクレジット表記に別段の記載がない限り、記事のクリエイティブ・コモンズ・ライセンスに含まれています。この記事のクリエイティブ・コモンズ・ライセンスに含まれていない素材で、あなたの意図する利用が法的規制によって許可されていない場合、あるいは許可された利用を超える場合は、著作権者から直接許可を得る必要があります。このライセンスのコピーを閲覧するには、http://creativecommons.org/licenses/by/4.0/。クリエイティブ・コモンズ・パブリック・ドメインの権利放棄（http://creativecommons.org/publicdomain/zero/1.0/）は、データへのクレジット表記に別段の記載がない限り、この記事で利用可能となったデータに適用されます。

転載と許可

この記事について
アップデートを確認する。CrossMarkで最新性と真正性を確認する。
この記事の引用
Liao、H., Ji、Y. & Sun、Y. ショートリードからの高分解能菌株レベルマイクロバイオーム組成解析. Microbiome 11, 183 (2023). https://doi.org/10.1186/s40168-023-01615-w

引用文献のダウンロード

受領
2022年8月24日

受理
2023年07月07日

出版
2023年8月17日

DOI
https://doi.org/10.1186/s40168-023-01615-w

この記事を共有する
以下のリンクをシェアすると、誰でもこのコンテンツを読むことができます：

共有可能なリンクを取得
コンテンツ共有イニシアチブSpringer Nature SharedItにより提供されています。

キーワード
菌株構成解析
メタゲノムデータ
k-mersインデックス構造
マイクロバイオーム
ISBN: 2049-2618

お問い合わせ
投稿に関するお問い合わせ: lyndie.manicani@springernature.com
一般的なお問い合わせ: info@biomedcentral.com
ブログを読む
BMCニュースレターを受け取る
記事アラートの管理
言語校正
著者のための科学的校正
ポリシー
アクセシビリティ
プレスセンター
サポートとお問い合わせ
フィードバックを残す
採用情報
BMCをフォローする
BMCツイッターページ
BMC Facebookページ
BMC微博ページ
このウェブサイトを使用することで、当社の利用規約、お客様の米国州プライバシー権、プライバシーステートメント、およびクッキーポリシーに同意したものとみなされます。プライバシーに関する選択/プリファレンスセンターで弊社が使用するCookieを管理する。

この記事が気に入ったらサポートをしてみませんか？