DNAシークエンスが開発されて以来、遺伝的変異に関する多くの研究が行われてきたが、エピジェネティックなレベルでの広範な研究は最近になって登場した。生体内のほとんどの細胞は、そのゲノム配列が同一であるが、組織や細胞の種類によっては、それぞれのアイデンティティを与えるエピジェネティックな修飾のパターンが異なる。DNAメチル化は最も重要なエピジェネティック・マークの一つで、主にCpGジヌクレオチドに生じる。ヒトゲノムには約2,800万個のメチル化部位が存在するため，450k個のアレイ（全CpGの1.6%しかカバーしていない）では，小さなdifferentially methylated regions （DMR）を検出するには不十分である[ref.1]。そのため、すべてのCpGメチル化レベルを適切に同定するには、データ量の多い全ゲノムバイサルファイトシーケンシング（WGBS）が必要となる。これらのデータセットを作成するためのコストは非常に高かったが、シーケンスコストの継続的かつ持続的な削減により、ますます多くのWGBSデータセットが作成されるようになり、包括的で再現性のある解析ツールが必要となっている。アラインメントやDMRの検出など、WGBS解析のさまざまな側面について、すでに多くのアルゴリズムが確立されている。しかし、適切なアルゴリズムを選択し、エンド・ツー・エンドの解析ワークフローに統合することは、可能なパイプライン設定がコンビナトリアル・エクスプロージョン（組合せ爆発）するため、簡単な作業ではない。エンドツーエンドのWGBS解析ワークフローを構築するには、入出力フォーマットや染色体の命名規則など、相互に作用するツールの要件が異なることがさらに大きな障害となる。これまでに開発されたエンド・ツー・エンドのパイプラインは、これらの問題をすでに考慮しており、ユーザーは生データと設定を提供するだけでよい。しかし、「結果と考察」の項で詳しく述べたように、これまでのアプローチでは、一般的な研究環境で必要とされる機能が欠けていたり、インストールの問題などの技術的な制限があったりした。その結果、これらの問題に対処するために、パイプラインアプローチを開発した。
ここでは、WGBSデータの自動インシリコ処理のためのワークフローであるwg-blimp（whole genome bisulfite sequencing methylation analysis pipeline）を紹介する。wg-blimpは、包括的なWGBSデータ解析パイプラインに加えて、データセットの検査を簡素化し、他の研究者と結果を共有する可能性のあるユーザーインターフェースで構成されている。論文図1は、提供される解析ステップの概要を示している。
FASTQファイルと参照ゲノムを入力として、wg-blimpはアライメントからDMR解析、セグメンテーション、アノテーションまでの完全なワークフローを実行する。アラインメントには、BWA-MEMの使用により、効率的でロバストなマッピングを提供するbwa-methを選択した。bwa-methはソフトクリッピングを使用して、一致しないリード部分をマスクするため、アライメント前のリードのトリミングは省略している。アラインメントからPicard toolkitを用いて重複を排除する。MethylDackelはbwa-methでの使用が推奨されているツールであるため、MethylDackelを用いてメチル化解析を行った。MethylDackelによって作成されたメチル化レポートに基づいて、wg-blimpはグローバルなメチル化統計を計算する。メチル化されていないラムダDNAは一般的にバイサルファイト処理の前にゲノムDNAに加えられるので、染色体ごとのメチル化の計算はオプションであり、C > T変換率の推定を可能にする。
品質管理（QC）には、FastQCを使用してリード品質スコアを評価する。全体および染色体ごとのカバレッジに関する情報を含むカバレッジレポートは，Qualimapによって生成される。Qualimapは、GCコンテンツ、重複率、クリッピングプロファイルなどのメトリクスもレポートするため、分析した各サンプルの詳細な品質評価が可能になる。Picard、Qualimap、FastQCによる品質レポートは、MultiQCを使用して1つのインタラクティブなHTMLレポートに集約される。
DMRコーリングには複数のアルゴリズムがサポートされており、metilene、bsseq、camelがよく使用されるツールである。複数のDMRコーリングツールを使用することを勧める。なぜなら、これらのツールは、重複するDMRのセットを識別するからである。ユーザーは、複数のパラメータを調整して、DMRコールの数を制御することができる。DMRコールの数を増やすと、通常、偽陽性のコールの割合が増えることになる。パラメータには、メチル化が異なると認識される領域に含まれるCpGサイトの最小数、比較する2つのグループ間の平均的なメチル化の最小絶対差、DMRあたりの最小カバレッジなどがある。metileneによるDMRコールには、Mann-Whitney U検定に基づくq値も含まれており、下流のフィルタリングに使用することができる。
さらに、非メチル化領域（UMR）と低メチル化領域（LMR）の検出を統合して、アクティブな制御領域を偏りなく特定する。このセグメンテーションはMethylSeekRを用いて実装されている。MethylSeekRは、ユーザーが定義した偽発見率（FDR）とメチル化カットオフのみを用いてモデルパラメータを自動的に推論する。MethylSeekRは，部分的にメチル化されたドメイン（PMD）と呼ばれる，非常に不規則なメチル化の領域の検出も行う。PMDの存在はUMR/LMRの検出に影響を与えており、しばしば事前には知られていない。その結果、wg-blimpはPMDの計算を伴う場合と伴わない場合のMethylSeekRワークフローを事前に実行する。MethylSeekRによって測定されたメトリクスに基づいて、ユーザーはUMRとLMRを分析するときにPMDを考慮するかどうかを決定することができる。
結果として得られたDMR、UMR、LMR、およびPMDは、EnsemblおよびUCSCデータベースで報告されているように、遺伝子、プロモーター、CpGアイランド（CGI）、および反復要素とのオーバーラップについてアノテーションされる。DMRごとの平均カバレッジはmosdepthを用いて計算し、カバレッジの低い領域のDMRコールをフィルタリングできるようにした。
wg-blimpパイプラインはワークフロー実行システムであるSnakemakeをベースにしている。Snakemakeは分析パイプラインのロバストでスケーラブルな実行を可能にし、障害が発生した場合に誤った結果を生成することを防ぐ。また，Snakemakeにはランタイムおよびメモリ使用量のロギング機能があるため，ボトルネックの探索やパフォーマンスの最適化が容易になる。ソフトウェアのバージョン変更に伴うエラーを最小限に抑えるため、依存関係の管理とインストールにはBiocondaを利用している。Biocondaの依存関係をAPIで変更するとパイプラインの機能が一時的に停止する可能性があるため、さらにwg-blimpのDockerコンテナを提供している。
分析ワークフローが完了すると、ユーザーは結果をwg-blimpのユーザーインターフェースに読み込むことができる。Rの機能をリアクティブなウェブアプリにシームレスに統合することができるR Shinyフレームワークを使用してインターフェースを実装した。このインターフェースは、QCレポート、パイプラインパラメータを集約し、呼び出し元の出力とアノテーションに基づいてDMRの検査とフィルタリングを可能にする。MethylSeekRによって計算されたUMRとLMRは、wg-blimpのShinyインターフェースを介してアクセスすることもでき、ユーザーはPMDを含めるかどうかを動的に選択することができる。ゲノムデータの可視化は、解析結果を確認する際によく用いられるため、IGV（Integrative Genomics Viewer）で使用するアライメントデータへのアクセスリンクも提供している。
インストール
#conda
mamba install -c bioconda -y wg-blimp
docker pull imimarw/wg-blimp:v0.9.7
> wg-blimp
# wg-blimp
Usage: wg-blimp [OPTIONS] COMMAND [ARGS]...
Options:
--help Show this message and exit.
Commands:
create-config Create a config YAML file for running the...
delete-all-output Remove all files generated by the pipeline.
run-shiny Start shiny GUI using configuration files for...
run-snakemake Run the Snakemake pipeline from command line.
run-snakemake-from-config Run the snakemake pipeline using a config file.
> wg-blimp create-config --help
# wg-blimp create-config --help
Usage: wg-blimp create-config [OPTIONS] FASTQ_DIR REFERENCE_FASTA GROUP1
GROUP2 OUTPUT_DIR TARGET_YAML
Create a config YAML file for running the Snakemake pipeline. Sample names
are either passed as comma seperated lists or are read from text files if
--use-sample-files parameter is set. Annotation files are automatically
downloaded if necessary.
Options:
--use-sample-files Load sample names from text files instead of
passing them as a comma-seperated list.
--genome_build [hg19|hg38|mmul10|None]
Build of the reference used for annotation.
--cores-per-job INTEGER The number of cores to use per job.
[required]
--help Show this message and exit.
実行方法
ここではdocker imageをランする。ラン後、分析用のshiny GUIにホスト側からアクセスするため、ポートをあけておく。docker内のwg-blimpコマンドは、デフォルトの設定では9898番を使う。ホスト側は9897番でhttp://localhost:9897にアクセスするなら、-p 9897:9898としてdocker imageをランする。ホスト側も9898なら-p 9899:9898。
docker run --rm -it -v $PWD:/data -p 9898:9898 imimarw/wg-blimp:v0.9.7
Githubには実際のデータのラン手順（link）が書かれている。wg-blimpは、デフォルトではIlluminaの命名規則でfastq名を判断する。wg-blimpをランするには、これらのfastqのパスとグループの関係を提供する。グループの定義は、fastqファイル名に基づいている必要がある；下の2群間比較のfastqファイルでは、blood1,blood2とsperm1,sperm2と指定している。
cd subsampled-blood-sperm/
wg-blimp run-snakemake fastq/ chr22.fasta blood1,blood2 sperm1,sperm2 results --cores=8
出力
alignment
methylation
segmentation
dmr
logs
qc
これらの結果はRのコンソールを使わずとも、R shinyでGUI環境で分析できる。出力ディレクトリにあるconfig.yamlを指定してwg-blimp run-shinyコマンドを実行する。
cd results/
wg-blimp run-shiny config.yaml
selectを選択。
statistics
Qualimap reportsとMultiQC reportsにもリンクしている。
parameters
DMRs
パラメータはリアルタイムで変更できる。
segmentation
他のコマンド
上記のようなwg-blimpを直接ランする”wg-blimp run-snakemake”コマンドの他、config.yamlを先に作成する”wg-blimp create-config”コマンドが用意されている。こちらは、configの修正が必要な時に使う。
wg-blimp create-config fastq/ chr22.fasta blood1,blood2 sperm1,sperm2 results config.yaml --cores-per-job 8
config.yamlが出力される。
configを開いてアノテーションのパスやパラメータを修正する（hg38以外は修正必須）。configのパラメータの詳細はGithub参照。
準備ができたらwg-blimpを”wg-blimp run-snakemake-from-config”コマンドでランする。。
wg-blimp run-snakemake-from-config config.yaml
引用
wg-blimp: an end-to-end analysis pipeline for whole genome bisulfite sequencing data
Marius Wöste, Elsa Leitão, Sandra Laurentino, Bernhard Horsthemke, Sven Rahmann & Christopher Schröder
BMC Bioinformatics volume 21, Article number: 169 (2020)
