教師ありディープラーニング法を用いたシス制御領域のゲノムワイド予測|BMC Bioinformatics

ディープラーニングは、アクティブエンハンサーとプロモーターをバックグラウンドから正確に区別する

ディープラーニングモデルが、エンハンサーとプロモーターを分離し、他の領域と活性状態間で区別する能力を調べました。我々は、手順を100回繰り返し、各よく特徴付けられた細胞型から私たちの（不均衡な）テストセットを予測するために私たちのバランスのとれたラベルトレーニングセット上の深いフィードフォワードニューラルネットワークを訓練しました。ディープモデルは、ゲノム領域上で実験的に導出された特徴を、確率を持つこれらの領域の入力および出力クラスラベルとして取ります（各クラスのサＡ−Ｅ、Ａ−Ｐ、Ａ−Ｘ、ｉ−Ｅ、ｉ−Ｐ、ｉ−Ｘ、およびＵＫとして、それぞれ、活性エンハンサー、活性プロモーター、活性エクソン、不活性エンハンサー、不活性プロモーター、不活性エクソン、およ活性Crrが転写を受けているという仮定の下で、活性は、CAGE転写開始イベントが焦点の組織で観察される領域に適用され、不活性は、他の組織で検出された領域を指すが、焦点組織では検出されない領域を指す。平均クラス単位の率を記録しました（すなわち、すべてのクラスの平均感度）、受信機動作特性曲線の下の面積（auROC）、および精度リコール曲線の下の面積（auPRC）。 1および追加ファイル1:図S1.

私たちが強調している結果には4つの側面があり、Crrのクラスと背景を区別するための教師付き深層学習アプローチの能力を確認しています。第一に、本発明者らは、活性増強剤と促進剤（Ａ−Ｅ対Ａ−Ｐ）とを区別することができる（図１０Ａ）。 1a）。我々は、それぞれ、正と負のトレーニングクラスとしてA-EとA-Pを使用しました。全体として、A-EとA-Pは非常に分離可能であることがわかりました。第二に、我々は、活性および不活性Crr（エンハンサーまたはプロモーターのいずれか）を区別することができます。図から。図1bおよび追加ファイル1：図S1Aでは、最大のトレーニングセットを有するGM12878、Helas3、Hepg2、およびK562上の平均auprcが0.95を超えており、エンハンサーおよこの論文の残りの部分では、限られたデータの可用性のために、ほとんどの分析でA549およびMCF7細胞株を除外します。第三に、予期せぬことではなく、不活性エンハンサーとプロモーターを区別することは困難である（追加ファイル1：図S1B）。八つの細胞型の平均クラスワイズ率の七つは0.80よりも低かった。不活性プロモーターの一部が存在するいくつかの機械を持っていることをいくつかの兆候があるが、それはそのような領域は、主に強力な転写因子結合第四に、我々は、I-E、I-P、A-X、I-X、およびUKをマージするスーパー背景（BG）クラスからA-EおよびA-Pを予測することの適用性をテストした（図4）。 1c）。六つの細胞型の結果は有望であった、すべて0.80auPRCを超えました。 A-EとA-Pがさらにマージされてスーパークラス（A-E+A-P）が形成されると、より高いパフォーマンスが達成されます（追加ファイル1：図S1C）。これらの六つの細胞型のすべてのauprcは0.89auPRCを超えていました。さらに、我々はまた、我々のラベル付けされたデータに、ランダムフォレスト法、別の最先端の分類器をテストしました。全ての実験設定で同様の性能が得られた。ランダムフォレスト法は、両方のエンハンサーの低い数を持っているA549とMCF7データセットのためにわずかに優れたパフォーマより多くの注釈付きエンハンサーが利用可能になることを期待して、我々はMLPを使用して、そのような畳み込みニューラルネットワークやリカレントニューラル

DECRESはFANTOM注釈領域に高い感度と精度を与えます

CRR予測のための私たちの教師付き深い方法の相対的な有用性を評価するために、我々は参照としてそれらは、真のゲノムの背景を反映した不均衡なセットで比較された。その結果を図１０で比較する。面積が大きく、より凸であるレーダーチャートを表示する2aは、より良い性能です。ラベル付けされた学習データが十分である場合、教師ありアプローチが好ましいことは直感的です。さらに、両方の監視されていない方法は、FANTOM5データの公開前に開発されたため、不利な立場にあります。しかし、これらの注釈はコミュニティによって広く使用されているため、標準に対するDECRESの相対的なパフォーマンスが重要です。全体として、DECRESはChromHMMと組み合わせた方法よりも優れており、同様の性能を提供しています。これらの監視されていない方法は一貫して活動的な増強物の検出のためのより低い感受性を有する(p=5。ＤＥＣＲＥＳ対Ｃｈｒｏｍｈｍｍのための５７Ｅ−５および９． 2b)および活性プロモーター検出のためのより低い精度(DECRES対ChromHMMおよびそれぞれ組み合わせた場合のp=7.36E-5および2.33E-4、two tailed Student’s t-test;図を参照。 2b）。 ChromHMMを使用すると、活性エンハンサーの感度は16.5％から48.4％の範囲であり（数値はエンコード予測エンハンサーのテストと一致しています）、deepモデルは69％（K562）から88.8％（GM12878）の範囲である。さらに、ChromHMMは49の最高の精密を達成する。活動的な促進者の予測のための8%、DECRESのための最高は84.3%であるが。

独立した実験データによるDECRES性能の評価

最初の評価は、教師ありモデルを訓練するために使用されるデータの種類であるCRRsのFANTOM eRNAベースの注釈に焦点を我々はさらにDECRESの性能を評価するために実験室で検証エンハンサーの二つの独立したコレクショ: Ｋ５６２細胞で試験された領域のＣＲＥ−ｓｅｑコレクションおよびＫ５６２細胞およびＨｅｐｇ２細胞で試験されたＭＰＲＡ（ｍａｓｓｉｖｅｌｙｐａｌｌｅｌｒｅｐｏｒｔｅｒａｓｓａｙ）コレクション。両方のインスタンスでは、発現を指示するために失敗した領域のセットが誤って評価された方法によって予測されるかもしれないが、また、実験手順データの性質を考えると、我々は実験ネガの一部が善意の調節領域であると予想している。

最初の独立したセットでは、予測されたK562エンハンサーと負の領域のサブセット（ChromHMMとSegwayの組み合わせ法によって予測される）は、CRE-seqを使用して実験室で評価された。その研究では、「結合された」予測された調節領域のわずか33％が実験で陽性であることが判明し、陰性セットの7％と比較した。従って、Ｋ５６２細胞の全ての利用可能な活性調節領域で訓練されたＤＥＣＲＥを使用して、２９８の対照領域と比較して、ＣＲＥ−ｓｅｑによって検証されたように、Ｋ５６２の活性エンハンサー活性を示す３８６の領域について、本発明者らの方法を検証した（追加ファイル１：表Ｓ３）。上記の結果と非常に一致して、実験的に検証された領域の65.5％（254/386）の感度が正常にA-Eとして予測された;残りの132領域が背景として予測された（どれも CRE-seq実験で不活性であった812試験された予測について、DECRESは53.3％（433/812）を陽性として分類した。２９８個の陰性対照領域について、ＤＥＣＲＥは、全てが陰性であると予測した（ＣＲＥ−ｓｅｑ実験において活性であった１６個を含む）。重要なのは、DECRESスコアが上昇するにつれて、予測の質が向上することです。我々は、DECRESによってA-Esとして予測された254と433実験的に陽性と陰性の組み合わせエンハンサーのDECRESメンバーシップスコアのヒストグラムを描いた（追加ファイル1：図S2）。分布は大きく異なります(p=0.014、両面マン-ホイットニーランクテスト）。

細胞選択的TFsの予測TF結合部位を含むK562とHepg2特異的な”強いエンハンサー”（ChromHMMによって予測される）は、超並列レポーターアッセイ（MPRA）を用いて試験された第二の独立エンハンサーの41％のみが有意に発現されることが検出された（p=0.05、両面マン-ホイットニーランクテスト）。ＭＰＲＡ陽性およびｍｐｒａ陰性エンハンサーのクラスを予測するためにＤＥＣＲＥＳを用いた。追加のファイル1：表S3の結果は、98.4％（120/122）と97を示しています。MPRA陽性エンハンサーの8％（182/186）は、それぞれK562とHepg2細胞のDECRESによってA-Esであると予測されたが、MPRA陰性エンハンサーの92.3％（179/194）と81.3％（217/267）はまだK562とHepg2のためのA-Esとして予測されたが、DECRESスコアの異なる分布（p=4.8E-6とp=2.3E-6k562およびhepg2それぞれ、両面マン-ホイットニーランクテスト）（追加ファイル1：図s2）。他の独立したデータと一致して、DECRESスコアが高いほど、それらが陽性である可能性が高くなります。

DECRESの性能に対するDNA配列特性の有用性の評価

最近の研究では、DNA配列特性がプロモーターおよびエンハンサーの認識、およびストリング配列カーネルこれは、プロモーター予測を改善するための機能としてCpG島を含めるための長い認識能力に基づいて構築されます。我々は、DNA配列の特徴は、プロモーターとエンハンサーの間、およびアクティブと不活性クラスの間を区別するために有益であり得るかどうかを決定私たちは、複数のシナリオで351個のシーケンス機能（当初は使用されていました）を使用してモデルを訓練しました。結果を図３に表示する。図3および追加ファイル1:図S3. まず、Ａ−ＥとＡ−Ｐを識別するための系列特徴に限定されたｄｅｅｐ法（図２）を用いて、ａ−ＥとＡ−Ｐを識別する。 3a)0.8567から0.9370までのAUPRCを送達し、配列属性が実際に有益であることを確認した。第二に、シーケンス機能は、論理的であるエンハンサーとプロモーターのアクティブと不活性状態を区別するための限られた有用性を持っている;実験的に導出された機能は非常にそれらを分離することができながら(p=1.90E-08と5.06E-08エンハンサーおよびプロモーターについては、両側スチューデントのt検定;図を参照してください。図3bおよび追加ファイル1:図S3A)。実験的特徴の非存在下で配列特徴を使用することは、８つの細胞型全てにわたってＡ−Ｅ、Ａ−ＰおよびＢＧを分類する際に低い性能を有する（ｐ＝１． 3c）。最後に、より良い結果は、実験的特徴と配列特徴を組み合わせることによって達成されなかった（図２のp＝2.79E-01、6.56E-01および1.17e-01）。 3、両側の学生のt検定）。

DECRESパフォーマンスのための主要な機能

実験データは時間がかかり、生成するために高価なことができるように、我々は計算の観点からCRR予測のための最も我々は、利用可能な72-135の機能を持っている四つの細胞型（GM12878、Helas3、Hepg2、およびK562）の二クラスと三クラス（A-E対A-P対BG）分類のためのランダム化された深

図4aおよび追加ファイル1:図S4Aは、三クラス分類のランダム化DFSおよびrandom forestによって検出された特徴重要度スコアを表示します。これらの方法によって生成される特徴重要度スコアは、異なる解釈をする必要があります。順方向選択と同様に、ランダム化DFSからの特徴重要度スコアは、スパースモデルの初期段階でどの特徴が優先されるかを反映し、ランダムフォレストによる特徴の重要度スコアは、他のすべての特徴との使用のコンテキストにおけるこの特徴の役割を示します。したがって、この研究では両方の方法を使用することで、データについて異なる洞察を得ることができます。我々の実験では、両方の方法は、すべての四つの細胞株にわたって重要性スコアによって示されるように、最も重要な特徴を捕捉することができ例えば、両方の方法は、Ｐｏｌ２、Ｈ３Ｋ４Ｍｅ１、Ｔａｆ１、およびＨ３Ｋ２７Ａｃが、ＧＭ１２８７８細胞株の背景から活性増強剤およびプロモーターを区別するのに有用であることにいくつかのケースでは、異なる尺度が互いに補完する。たとえば、h3K4Me2とH4K20Me1は、ランダム化DFSによって主要な機能としてマークされており、追加ファイル1の箱ひげ図で示されているように説得力があります：図S4Bと図S6-S13ですが、random forestによって見落とされています。 Tbpは、GM12878およびHelas3細胞におけるランダムフォレストによって強調表示されたが、無作為化DFSによってピックアップされませんでした。追加ファイル1:図S6およびS7でこの機能の箱ひげ図を調べると、この機能はアクティブエンハンサーとプロモーターを背景から区別するために差別的であることが明らかになりましたが、アクティブエンハンサーとプロモーターの間には劇的な違いはありません。ランダムフォレストモデルに組み込まれた重要なフィーチャは、DFSプロセスの後の段階まで組み込まれない場合があります。例えば、Ｋ５６２細胞株において、Ｃ−Ｍｙｃは、Ｒａｎｄｏｍｆｏｒｅｓｔによって強調され、これは、追加のファイル１：図Ｓ１２に示されるように、実際には妥当であり、ＤＦＳプロセスの初期

ゲノム注釈における機械学習法の開発のために、必要な特徴の数を最小限に抑えることは、コストを削減し、生物学的解釈のための能力を増加させる。図4bおよび追加ファイル1:図S5Bは、選択された特徴の数がそれぞれ3クラスと2クラスの分類で増加するにつれて、テストauprcの変更を示していどちらの場合も、テストauPRCsは最初の特徴、そして性能のプラトーのために劇的に増加します。ランダム化されたDFS曲線とランダムフォレスト曲線を比較すると、単一の最適曲線が存在しないことがわかります。いくつかの主要な機能は、良好な予測性能のために十分である。必要とされる最適な数の特徴を定義するために、図1の曲線を近似します。図4bおよび追加ファイル1:図S5Bで、ランダム化されたDFS曲線上の傾きが0.5の線の交点を選択しました(方法を参照)。 A-E、A-P、および背景（10の特徴）を区別することを意図した三クラスモデルと比較して、二クラスCRR予測（6つの特徴）に必要な特徴が少ない。

3つのクラス予測(A-E,A-P,BG)の上位10個の特徴の分布は、追加ファイル1:図S4Bに示されています。各セルの上位10個の特徴を使用して、GM12878,Helas3,Hepg2,K562でそれぞれ0.9022,0.9156,0.8651,0.8565のauprcが達成されました。これらのトップ機能の半分は、h3K4Me1、H3K4Me2、H3K4Me3、およびH3K27Me3は、既存の知識と一致して、一般的に三クラスモデルのための機能を選択したヒストンの修正です。転写因子（共因子を含む）の中で、TAF1およびp300、ならびにRNAポリメラーゼII（Pol2）が頻繁に選択され、これも既存の知識と一致している。

追加ファイル1:図S5Cは、選択された上位6つのフィーチャの箱ひげ図を、2つのクラス予測のランダム化DFSによって示しています。これらの機能を使用して、0.9561、0.9627、0.926、および0.9555のauprcは、それぞれ、四つの細胞型で得られました。ほとんどのフィーチャでは、a-EおよびA-Pの値の範囲は、背景カテゴリに対して高くなります。選択された特徴の半分は、H3K4Me2、H3K27Ac、およびH3K27Me3を含むDNase-seqおよびヒストン修飾チップ-seqデータである。これらの特徴の箱ひげ図は、A-Eとa-Pを背景から区別することを示しています。

DECRESのゲノム全体の予測の大部分は、他の方法でサポートされています

我々は、6つの細胞型（A549とMCF7は除外された）のゲノム全体のCrrを予測するために、すべての参照（ラベル付けされた）データを使用して、2クラスと3クラスの多層パーセプトロン（MLP）モデル（メソッドを参照）を訓練しました。 2クラスモデルは、ゲノムの4.8％を占める227,332Crr（隣接領域がマージされた）を同定した（追加ファイル1：表S4）。 9153Crrの合計は、すべての六つの細胞型にわたって遍在的に予測されました。 3クラス予測のために、我々は301,650A-E領域（ゲノムの6.8％）と26,555a-P領域（ゲノムの0.6％）と一緒に11,886ユビキタスA-Esと3678ユビキタスA-Psを得た。 6つの細胞型すべてのゲノムワイド予測は、追加のファイル2で利用可能です。

次に、GM12878、Helas3、およびK562の組み合わせ予測およびdReg予測と予測Crrの重複を調べました。ＤＥＣＲＥＳによって予測されるＣｒｒの大部分は、Ｇｍ１２８７８、Ｈｅｌａｓ３、およびＫ５６２について、それぞれ、結合またはＤＲＥＧのいずれかからの結果、具体的には８６.１３％、７６.１３％、および８３.６３％と重 5). DECRES予測のサブセット(GM12878では13.87%、Helas3では23.87%、K562では16.37%)は、他の二つのツールからの予測と重複しません。特に、結合された予測の大部分（Helas3の56.78％、GM12878の55.99％、および36。36%on K562)は、監視された方法からのものと重複しません,その低い観察された検証率と一致しています. さらに、DECRES予測は、A-P領域とA-E領域の両方でより細かい分解能を持つ傾向があります（例については、追加のファイル1：図S14を参照）。

私たちは、VISTAエンハンサーセットによってサポートされている私たちのゲノム全体の予測の中でどのように多くの調査しました。 VISTAのエンハンサーの大半は非常に開発全体で保存されているという事実にもかかわらず、我々はまだ実験的に確認され、未確認のVISTAのエンハンサーの37.1％（850/2,293）は、予測されたA-Esと重複し、これらのVISTAのエンハンサーのわずか4.8％（110/2,293）は、予測されたA-Psと重複していることがわかります。実験的に確認されたVISTAエンハンサーの結果は似ています（482/1、196=40.30％と60/1、196=5.02％がそれぞれA-EsとA-Psに重複しています）、これは私たちの予測されたアク私たちの予測と重複していないVISTAエンハンサーの割合は、開発中または私たちのフォーカス細胞株以外の他の細胞型で特異的に活性である可能性があ

DECRES extends THE FANTOM enhancer atlas

ernaのケージ信号の深さが限られているため、活性（または転写された）エンハンサーの一部はenhancer atlasの元のコンパイルでは検出されません。したがって、我々はeRNA信号が元のアトラス閾値設定を下回っていたため、追加の部分的にサポートされているエンハンサーを識別しようとした。前回の研究では、200,171双方向転写（BDT）遺伝子座の合計は、808細胞型および組織のケージタグを使用して、ヒトゲノム全体で検出されました。エクソン内のBDT遺伝子座を除外した後、部分的にサポートされている102,021のBDT領域のセットが残り、そのうち43,011のバランスの取れた遺伝子座（両側の同様のeRNAレベル）がFANTOM enhancer atlasを構成する。よりアクティブなエンハンサーの候補は、六つの細胞型のそれぞれのために検出することができるかどうかを調べるために、我々はそのアクティブなアトラス領域上のMLPを訓練し、すべての102,021BDTサイトのためのクラスを予測しました。 102,021個のBDT遺伝子座のうち、ほとんどが特定のセル内の負の領域として分類された（追加ファイル1: 一方、平均して１３，３１６がＡ−Ｅｓとして予測され、８３４のみが細胞型当たりＡ−Ｐｓとして予測された。元のエンハンサーアトラスにおける不活性なエンハンサーのかなりの数（平均6535）は、BDTデータが任意の所与のサンプルのために不完全であるという仮定と一致して、私たちのモデル（追加ファイル1：表S6）によって活性として予測された。元のアトラスによって除外された平均５５１４個のＢＤＴ遺伝子座を、細胞型ごとのＡ−Ｅｓとして予測した。解析された六つの細胞型の上に、38,601BDT遺伝子座の合計は、16,988は、元のFANTOMエンハンサーアトラスの拡張を表す-Es（追加ファイル3）として予測されました。元のFANTOMエンハンサー atlasから21,398 43,011エンハンサーのうち、ここで分析された六つの細胞で活性として予測されていないが、これらの領域は、分析するのに不十分な特徴がある他の802細胞で活性である可能性があることに注意してください。

機能的およびモチーフ濃縮解析を用いたDECRES予測の計算検証

GREATを用いてゲノム全体の予測A-EsおよびA-Psに対して機能的濃縮解析を行った。 GM12878細胞の場合、予測されたエンハンサー領域の79％が遺伝子Tssから5キロ塩基対（kbps）以上離れており（追加ファイル1：図S15A）、予測されたプロモーターの47％が注釈された遺伝子Tssに対して5kbps未満である（追加ファイル1：図S15B）。残りの５つの細胞型についても同様の統計が得られた。 GM12878特異的Crrの注釈分析は、近位遺伝子がに関連していることを示している：遺伝子オントロジー（GO）注釈（追加ファイル1：図S15C）からの免疫応答;MSigDB経路注釈（追加ファ結果は、細胞のリンパ芽球様系統と一致する。次に、FANTOM enhancer atlas（「not in atlas」）で以前に報告されていないBDT支持予測エンハンサーについて機能濃縮分析を行った。結果は上記の分析と完全に一致しています（追加ファイル1：図S16）。

我々はさらにHOMERを用いて予測された細胞特異的Crrとnot-in-atlasエンハンサーのモチーフ濃縮解析を行った。予測された領域は、JASPAR結合プロファイル（追加ファイル1：図S15Fと図S16-S26）の両方の一般的な細胞プロセスを維持するTFsと細胞関連機能における選択的な役割を持つTFsに関連付けられているようなモチーフのために濃縮されている。例えば、Jun-、Fos-、およびEts関連因子のモチーフは、すべての六つの細胞型からの領域に富化された。これらのTFsは、分化、増殖、またはアポトーシスなどの一般的な細胞の進行を調節する。細胞に適したTF濃縮物を各細胞について観察した（追加のファイル１：表Ｓ７に要約）。例えば、造血において重要な役割を果たすRUNX1および他のRunt関連因子は、GM12878で観察される（追加ファイル1：図S15Fおよび図S16）。免疫応答および炎症応答に関与する遺伝子を調節するC/EBP関連因子は、子宮頸部で発現される（追加ファイル1：図S17およびS18）。 HNF1A、HNF1B、FOXA1、FOXA2、HNF4A、およびHNF4G因子は、肝臓特異的遺伝子を調節する（追加ファイル1：図S19およびS20）。ＮＦＹ因子は、ｇａｔａ１と協力して、ｋ５６２における赤血球特異的転写を媒介する（追加ファイル１：図Ｓ２５およびＳ２６）。

組み合わせた方法からA-EおよびA-P予測について機能分析および濃縮分析を行い、その結果を追加ファイル1：図S27-S30に報告しました。組み合わせた方法によって予測されるプロモーターの大部分は、エンハンサーに類似している既知の遺伝子TSSsの遠位である。例えば、細胞株ＧＭ１２８７８上では、ｄｅｃｒｅｓプロモーターの４７％と比較して、結合されたプロモーターの２２％のみが、注釈された遺伝子Ｔｓｓに対して５ｋｂｐ未満に位置する。さらに，複合法によって予測されたＣｒｒｓの機能解析は，ＤＥＣＲＥＳ予測よりもＧＯ生物学的プロセス，Ｍｓｉｇｄｂ経路，および疾患オントロジーの有意項がはるかに少ないかゼロであった。両方の方法のモチーフ解析結果は一貫しています。

教師ありディープラーニング法を用いたシス制御領域のゲノムワイド予測