GREINの概念図を図1に示す。 個々のRNA-seqデータセットはGREP2パイプラインで処理され、R Expression Setsとしてローカルに保存される。 ユーザーはGREINのグラフィカルユーザーインターフェイス（GUI）から前処理済みのデータセットにアクセスして解析したり、まだ処理されていないデータセットを処理に提出したりすることができる。 GUIによるワークフローは、データの解析と可視化、統計解析、転写シグネチャーの構築、および差次的発現遺伝子のシステムバイオロジー的解釈を容易にします。 GREINとバックエンドパイプライン（GREP2）は共にRで書かれており、それぞれDockerコンテナとRパッケージとしてリリースされています。 GREINのグラフィカルユーザーインターフェイスは、Rで動的なWebアプリケーションを構築するためのWebフレームワークであるSiny16で実装されており、https://shiny.ilincs.org/greinのWebインスタンスは、負荷分散されたSinyサーバーの堅牢なDocker群によって展開されている。
GREP2の概略ワークフロー、Webインターフェース、GREINの出力。 GEOデータセットはGREP2パイプラインで系統的に処理され、バックエンドのデータセットライブラリに格納されます。 GREINのGUIワークフローは、処理されたデータセットの包括的な解析と可視化を容易にします。
GREINのユーザーフレンドリーなGUIワークフローは、品質管理の検討やデータセット全体における発現パターンの可視化、将来の研究の実験デザインに関する情報を得るためのサンプルサイズやパワー解析、統計的遺伝子差発現、遺伝子リストのエンリッチメント、ネットワーク解析などのRNA-secqデータの代表的な再利用シナリオを促進します。 遺伝子発現差解析モジュールは、標準的な2群比較に加えて、共変量やバッチ効果を考慮した一般化線形モデルの適合もサポートしています。 インタラクティブな可視化・探索ツールとして、クラスター分析、インタラクティブヒートマップ、主成分分析（PCA）、t-distributed stochastic neighbor embedding (t-SNE)などが実装されている。 (補足表S1). また、MetaSRAプロジェクト13が提供するヒトRNA-seqサンプルやデータセットのオントロジーアノテーションを検索することもできる。 処理された各ヒトRNA-seqサンプルには、Disease Ontology、Cell Ontology、Experimental Factor Ontology、Cellosaurus、Uberonなどの生物医学オントロジーのMetaSRAマッピングが付されています。 遺伝子発現の差の生物学的解釈は、遺伝子リストやパスウェイの濃縮分析、差次的発現（DE）遺伝子のネットワーク分析などの典型的なポストホック分析を実行する他のオンラインツールへの直接リンクによって支援される。 これらの分析ウェブサービスへの接続は、遺伝子発現の差分シグネチャー（すなわち、分析された全/上/下制御遺伝子に対する遺伝子発現の平均変化と関連するp値のリスト）をiLINCS17 (Integrative LINCS) に送信することで実装される。iLINCSは、最近リリースされた Connectivity Map L1000 signatures18 のシグネチャー結合性分析も提供している。 GREINの解析ワークフローの詳細については、補足資料およびGREINの「ヘルプ」を参照されたい。
Key functionalities
Search or submit for processing
ユーザーは「Search for GEO series (GSE) accession」ボックスで既に処理済みのGEOデータセットを検索し、処理されていないデータセットについては処理用にデータを提出できる（補足図 S2）。 現時点では、GEOのヒト、マウス、ラットRNA-seqデータセットの大部分は前処理済みであり、GEOデータセットのユーザ提出はごくまれにしか必要ないだろう。 ユーザーは’Processing console’タブでリクエストしたデータセットの処理状況を確認することができる（Supplementary Fig.） その他、データセットのメタデータによるキーワード検索や、MetaSRAオントロジーアノテーションによる生物医学オントロジー検索が可能である。 GREINには、発現パターンを可視化するためのインタラクティブでカスタマイズ可能なツールが付属している。例えば、クラスター化した遺伝子やサンプルのインタラクティブなヒートマップ、全サンプルまたはサブセットの密度プロット、PCAやt-SNEなどの2次元および3次元次元解析と可視化によるグループ間およびグループ内の変動解析などがある（図2）。 5673>
GREINにおける探索的分析プロット。 (A)相関ヒートマップでは、細胞株内の相関が高く、細胞株間の相関が低いことがわかる。 一般的に各細胞株内での相関が高いことは、転写プロファイルの質が高いことを示しています。 (B) 変動性の指標として中央値絶対偏差を用い、最も変動性の高い遺伝子上位500個のピアソン相関に基づく階層的クラスタリング。 データは正規化され、平均値でセンタリングされている。 (C)細胞株の3次元主成分分析プロット。 (D) 処理条件と細胞株の2次元t-SNEプロットは、細胞株の明確な分離を示し、次にRNA分画は、RNA-seqプロファイル間の変動の2つの支配的な原因を示す。 Deelenらによる最近の研究19では、処理された65,000のパブリックRNA-seqサンプルの半分以上が、QC問題のために削除されなければならなかった。 GREINはサンプルを削除するのではなく、各サンプルの生のシーケンスデータとシーケンスマッピングの包括的な品質管理（QC）レポートを提供し（補足図S7）、ユーザーはどのサンプルをダウンストリーム解析から除外すべきかを決定することができる。 RNA-seqデータの再解析において、類似の生物学的サンプルを用いた将来の研究に対する適切なサンプルサイズの見積もりは、しばしば重要な動機付けとなることがあります。 検出力分析はまた、現在のデータセットにおける偽陰性率のポストホック分析も容易にします。 統計的検出力の不足や遺伝子間の統計的検出力の差は、誤った結論につながる偽陰性を生み出す可能性があります20。 Power curve “セグメントは、単一遺伝子に基づく異なるサンプル数に対する検出力の推定値を提供する（Fig. 3A）。 ユーザーはパラメーターのデフォルト値を変更することができる。 Detectability of genes’ plotは、選択されたグループと遺伝子ごとの分散に基づいた各遺伝子の検出力推定値を表示する(Fig. 3B)。 5673>
Power analysis for assessing transcriptional changes in non-malignant MCF10A cell line.（「非悪性腫瘍細胞株における転写変化の評価」）。 (A)各グループのサンプル数が異なる場合の単一遺伝子ベースの検出力推定値は、最小倍率変化を2、統計的有意性をα=0.01とした。 (B)FDR≦0.1、各群2サンプルでのlog2CPM-BCOV平面上での遺伝子単位の検出力
Differential gene expression
RNA-seq 実験では遺伝子発現シグネチャの差の作成と解釈は典型的な分析場面である。 GREINでは、実験共変量やバッチ効果を調整した、あるいは調整しない2つのサンプル群間の遺伝子発現を比較することにより、シグネチャーを作成することができます。 GREINは、グループやサブグループの再配置、特定のサンプルの選択など、柔軟な対応により複雑な実験デザインにも対応することができる。 発現シグネチャーの可視化には、FDR（False Discovery Rate）順位による発現量上位遺伝子のヒートマップ（補足図S15）、log fold change vs. log average expression (MA) plot（補足図S16）、遺伝子検出性 plot（補足図S17）などのインタラクティブなグラフィックが利用可能である。 また、Differential expression signatureは、false negativeの可能性がある結果を考慮してもしなくても、濃縮と連結性解析のためにiLINCSに直接エクスポートすることができる。 非悪性乳房上皮およびトリプルネガティブ乳がん細胞株における低酸素の転写および翻訳制御の解析
最近公開されたGEO RNA-seq データ（GSE104193）の再解析により、GREINの使用法を実証します。 Seséら21は、低酸素とmTOR（mechanistic target of rapamycin）阻害剤治療の組み合わせでホルモン不応性トリプルネガティブ乳がん（TNBC）サブタイプの転写および翻訳制御を検討した。 特に、TNBC（MDA-MB-231）と非悪性乳上皮（MCF10A）細胞を正常酸素（21％O2）と低酸素（0.5％O2）条件下に曝し、および/またはmTORC1および-2阻害剤PP242で処理して発現プロファイルを分析した。 各サンプルについて、total (T) および polysome-bound (P) mRNA の配列が決定された。 5673>
GREIN で処理したデータセットを探索的に解析した結果 (図 2)、サンプル間の最も大きな変動要因は、2 つの細胞株の違いに起因することがわかりました。 このことは、発現プロファイルの相関解析（図2A）、絶対値偏差の中央値に基づく上位500個の変動遺伝子の階層的クラスタリング（図2B）、サンプルの3次元PCAプロット（図2C）、2次元t-SNEプロット（図2D）により再強化されている。 さらに、同一細胞株における発現プロファイル間の高い相関（図2A）は、遺伝子発現測定における良好なS/N比を示す。 2次元t-SNEプロットで示されたデータの追加的な部分構造は、異なる属性に従ってサンプルをペイントすることで調べられた（補足図S1）。 この解析により、各細胞株内の分離は、異なるmRNA分画によって誘導され、次に実験条件間の差によって誘導されることが明らかになった。
次に、このデータセットで観察された生物学的変動のパターンに基づいて、GREINを用いて統計力解析を実施した。 低酸素にさらされた各細胞株の転写プロファイルを、PP242で処理するかしないかで検討し、4つの比較を行いました。 グループ間で少なくとも2倍の発現差があると仮定し、α=0.01の統計的有意性で、各グループに2つの複製しかない場合、差次的発現として検出される遺伝子の統計的検出力はすべての比較で0.55以下である（表2）。 その結果、2倍の発現変化を検出するために80%の検出力を得るには、1群あたり4回の複製が必要であることがわかった（表2、図3A）。 典型的なRNA-seq実験では、ほぼすべての遺伝子について遺伝子発現を定量化するには、2000万から3000万の配列深度で十分である4,22が、このデータセットでも明らかである。 また、各遺伝子が差次的に発現していると検出される統計的な力を「遺伝子の検出力」プロットから評価した。 遺伝子の平均log of counts per million (CPM) 値を遺伝子ごとの生物学的変動係数 (BCOV) に対してプロットし、対応する遺伝子について検出力を計算した (Fig. 3B)。 統計的有意性を推定するために、0.05の制御された偽発見率と10%の期待される真陽性の割合を使用した。 検出力が0.8以上の場合、その遺伝子は低酸素状態で差次的に発現していると検出できると定義した。 予想通り、BCOVとpowerの間には逆相関が存在する（Fig.3B）。 また、遺伝子の差次的発現を検出するための検出力は、log CPMまたは効果量が大きいほど増加する。
研究の目的の一つは、MCF10AおよびMDA-MB-231細胞株の両方でPP242処理あり、なしの低酸素および標準酸素条件における転写の変更を分析することであった。 各細胞株について、「複製」を共変量として扱うことでバッチ効果を調整しながら、コントロールサンプルに対する低酸素および低酸素＋PP242サンプルそれぞれの差分発現解析を行い、総mRNAにおける低酸素および低酸素＋PP242の転写シグネチャーを作成しました。 低酸素および低酸素＋PP242のいずれにおいても、MCF10A細胞株ではMDA-MB-231と比較してより多くの遺伝子に差次的発現（DE）が認められた（図4A）ことから、おそらく腫瘍細胞株は低酸素に対処する能力がより優れていることが示唆された。 また、この解析では、ほとんどの非差異的発現遺伝子も検出されないことから、偽陰性結果である可能性があることが示された。 これは、平均的なBCOVで差次的発現遺伝子を一貫して同定するためには、1群あたり4サンプルが必要であることを示した検出力解析と一致する。 5673>
遺伝子の発現差と検出可能性。 (A）総mRNA分率について正常酸素との比較で、差次的に発現しておらず検出できない（NDE&NDT）、差次的に発現している（DE）、差次的に発現していないが検出できる（NDE&DT）遺伝子数（log10スケール）を示した。 検出力≧0.8の遺伝子を検出可能（DT）、FDR < 0.05の遺伝子を分化可能（DT）と呼ぶことにした。 (B) 最初の比較（MCF10Aと低酸素）の遺伝子検出性プロットで、上記の遺伝子リストをそれぞれのフォールド変化（FC）と共に可視化した。
影響を受ける生物学的経路という観点から発現差分を解釈するために、iLINCSを介してオンラインの濃縮ツール（DAVID23、ToppGene24、Enrichr25、Reactome26）に低酸素の遺伝子発現差分のサインを提出した。 提出されたシグネチャーには、真陽性と真陰性の可能性が高いDEとNDE&DT遺伝子の複合リストが含まれている。 遺伝子は、統計的検出力とFDRのそれぞれ0.7と0.01のカットオフ値に基づいて選択された。 図5は、MCF10低酸素シグネチャーのToppGeneから得られた濃縮結果を示している。 有意に濃縮された（FDR < 0.05）ToppGeneおよびDAVID機能アノテーションツールによる遺伝子オントロジー（GO）カテゴリー上位10位には、低酸素への応答、酸素濃度低下への応答、血管新生、細胞増殖の制御、酸化還元過程、および生物学的刺激への応答が含まれ、両細胞株で共通している（補表S2および補表S3）。 これらのカテゴリのほとんどは、オリジナルの研究と一致している。 さらに、ToppGeneスイートは、両方の細胞株で活性化された低酸素誘導因子（HIF-1-α）転写因子ネットワークを特定した（補足表S4および補足表S5）。
iLINCSを介してToppGeneからの重要な経路および遺伝子オントロジー（GO）カテゴリーのいくつかのスナップショットを示す。 これらのカテゴリは、DEとNDE&DT遺伝子の複合リストを使用して、MCF10A細胞株の低酸素と正常酸素の比較で発見されたものである。 5673>
最後に、iLINCSとのGREIN接続を利用し、アップロードしたシグネチャーをLINCS27コンセンサス（CGS）遺伝子ノックダウンシグネチャーと「接続」させた18。 その結果、アップロードしたシグネチャーと有意に（p値< 0.05）接続したLINCSコンセンサス遺伝子ノックダウンシグネチャーが3,727件見つかった。 さらに、上位100のシグネチャーの標的遺伝子を選択して、エンリッチメント解析を行った。 その結果、低酸素に対する細胞応答と酸素によるHypoxia-inducible Factor（HIF）の制御が、両細胞株の活性化上位10経路のリストに含まれていることがわかった（補足表S6および補足表S7）。 この解析では、最初の濃縮解析と同様の機能カテゴリーが得られたが、検出力解析によれば十分に高発現であるにもかかわらず差次的に発現しないいくつかの標的遺伝子を示唆することによって、最初の解析を補完するものであった。 これらの2つの結果を結びつけると、これらの遺伝子は低酸素に対する応答のより高いレベルの制御因子である可能性が示唆される。