An easy-to-use, autorun package for DNA methylation analyses
DNAメチル化データ解析をより便利に完了するために、すべての機能をパッケージ化してDNAメチル化解析用の簡単で自動実行できるパッケージとしました。 図1にBatMeth2の主な特徴を示す。 1) BatMeth2は、効率的で正確なアライメント性能を有しています。 2) BatMeth2は、個々のシトシン部位や、染色体全体、遺伝子領域、トランスポーザブルエレメント（TE）など、任意の機能領域のDNAメチル化レベル（ML）を算出することができます。 3) 異なる統計アルゴリズムを統合した後、BatMeth2は、任意の領域、任意の数の入力サンプルおよびユーザーの要求に対して、差分DNAメチル化解析を実行することができます。 4) BS-Seq データの可視化（染色体上および遺伝子上の DNA メチル化分布）と差分メチル化アノテーションを統合することで、BatMeth2 は DNA メチル化データをより明確に可視化することができます。 BatMeth2ツールの実行中に、サンプルの統計情報のためのhtmlレポートが生成されます。 サンプルhtmlレポートの詳細をhttp://htmlpreview.github.io/?https://github.com/GuoliangLi-HZAU/BatMeth2/blob/master/BatMeth2-Report/batmeth2.html.
BatMeth2のワークフローに示している。 2つの大きな矢印は入力ファイルまたは出力ファイルを意味する
BatMeth2 has better mapping performance on simulated BS-Seq data
我々はまず、75塩基対（bp）、100bp、150bpからなるリードと異なるバイサルファイト変換率（0から100％まで、ステップ10%）のシミュレーションデータセット（インデルを含まない）ですべてのアライナーを評価しました。 これらのデータセットは、FASTX-mutate-tools、wgsim (v0.3.0) およびSAMtools (v1.1) のシミュレータを使用して、ヒトゲノム (UCSC hg19) からシミュレートしました。 シミュレートしたリードを参照ゲノムにマッピングし、最大で2つのミスマッチを許容した。 シミュレーションしたリードの元の位置が分かっているので、マッピング出力と元の位置を比較することで、すべてのプログラムの精度を評価することができました。
FASTXとwgsimからの異なるリード長のシミュレーションデータセットによるBS-Seqアライナ全製品の評価。 ビスルファイト変換率の異なるシミュレーションデータは、異なる形状で表示されています。 異なるアライナーの結果は、シンボルの色を変えて表示しています。 各パネルの左上付近の結果は、正しくマッピングされたリードの数が多く、誤ってマッピングされたリードの数が少ないことを表しています。 弊社のアライナーBatMeth2の結果は、異なるシミュレーションBisulfiteデータセット
において、ミスマッチとインデルが混在する一般のシミュレーションBSリードをアラインした場合、wgsim-シミュレーションインデルアベラントデータセットの結果は、他の手法に比べてBatMeth2が良い性能（2番手アライナより1~2%良い）を示していることが分かります。 BS変換率の上昇に伴い、全てのソフトのアライメント精度が低下していることが分かります。 このような異なる条件下では、BatMeth2がより良いパフォーマンスを発揮しています。
BatMeth2 has better mapping performance on real BS-Seq data
BatMeth2 の性能を実際の BS-Seq データセットで試すために、ペアエンド BS-Seq データセットをダウンロードして、SRA SRR847318 から2 × 90 bpペアエンドリード100万本をランダムに抽出しました。 SRA SRR1035722から100万本の2×101 bpペアエンドリード、SRA SRR3503136から100万本の2×125 bpペアエンドリードを評価用として使用しました。 これらのデータセットは健康な細胞株や組織から得られたものであるため、構造変異の数は少ないと予想されます。 そこで、ペアエンドのシングルエンドリードを用いてこれらの実データをアライメントし、ペアアライメントから得られたマッピングの一致率と不一致率を評価し、アライメントの正解率と不正解率を見積もりました。 ペアエンドリードの挿入サイズは約500bpであるため、一対のパートナーリードが公称距離500bp以内にマッピングされればコンコードとみなされ、そうでなければ、一対のパートナーリードはディスコードとみなされる可能性があります。 図3
異なるアライナによるペアエンドリードでのアラインメントの一致率および不一致率を示している。 実際のバイサルファイトシーケンスリードについて、マッピング品質が高いものから低いものまで、一致するアラインメントと不一致するアラインメントの累積カウントを表示します。 BSmapとbowtieベースのアライナーは、マッピング品質のスコアがないため、別々に1点のみ存在します。 Bismark-bowtie2L15は、種長15
のbowtie2アライメントを意味する。さらに、表1に各プログラムの相対実行時間を示した。 デフォルト設定のBatMeth2は、ほとんどの公開アライナーよりも高速に動作し、BWA-methとBatMethに匹敵するものであった。 Bismark2 (Bowtie2を基本マッピング法とする)、BS Seeker2、BSmoothはより長い実行時間を必要とします。
DNA methylation calling
DNA methylation callingのソフトウェア間の精度を評価するために、ENCODE (Encyclopedia of DNA Elements) からIMR90セルラインの450K bead chip dataをダウンロードした。 また、ENCODEからIMR90細胞株の全ゲノムバイサルファイトシーケンス（WGBS-Seq）データ（42.6Gbases）をダウンロードした。 各ソフトウェアについて、WGBS-Seqのリードをアライメントし、DNAメチル化レベルを算出した。 その後、450K Bead Chipデータの同一部位のMLと比較した。 ソフトウェアによるWGBS-SeqデータからのDNA MLと450 K Bead ChipからのMLの差が0.2未満の場合、コール結果は正しいと定義し、それ以外は不正解とした
その結果を表2に示す。 全ソフトウェアの正しい結果の重複を追加ファイル1に示す。 図 S2 に示す。 BatMeth2とBiscuitは同程度の性能を持ち、他のソフトよりも優れていることが分かる。 結論として、BatMeth2はBSリードのアライメントとDNA MLコールの両方の精度を向上させる。
BatMeth2 は可変長インデルを許容しながらBSリードをアライメントする
がんには健康細胞より著しく高い割合でインデルが含まれています。 そこで、BatMeth2が長さの異なるインデルを含むBSリードをアライメントできるかどうかを検証するために、ENCODEからHepG2（肝臓肝細胞がん、がん細胞株）のWGBSデータ（75Gbases）と450K Bead Chipデータをダウンロードした。 HepG2 WGBS-Seqデータのアライメントを行った後、リードにおけるインデル長分布を確認した。 追加ファイル1: Figure S3Aによると、検出されたインデルの長さは主に1 bp~ 5 bpの範囲に分布し、最も長いインデルの長さは40 bpであった。 また、アラインメントリードの2.3%にインデルが含まれていることが確認された。 この結果から、BatMeth2は様々な長さのインデルを含むリードをアライメントできることが分かった。
次に、DNAメチル化コールに対するインデル検出の効果を検証した。 BatMeth2については、HepG2のデータに対して、インデル検出あり、なしの2つのオプションを実行しました（つまり、BatMeth2の-Iパラメータを設定）。 また、Bismarkにはインデルコール機能がないため、インデル検出ありのDNAメチル化コールのリファレンスとして、HepG2のWGBS-Seqデータに対してBismarkを実行しました。 BatMeth2とBismarkのDNAメチル化コールと450K Bead Chipデータからのコールを比較した。 その結果をAdditional file 1に示す。 BatMeth2-noIndel “は、インデルを検出しないBatMeth2に相当します。 インデル検出を行わない場合、BatMeth2の結果はBismarkの結果（Bowtie1を基本マッピング手法とした場合）よりもわずかに優れていることが分かります。 しかし、BatMeth2の結果は、Bismark（Bowtie1を基本マッピング法とする）よりもわずかに優れているに過ぎません。 さらに、BatMeth2はBatMeth2-noIndelやBismark（Bowtie1）よりも多くのDNAメチル化部位を検出できることが分かります。 なぜBatMeth2の方がインデル検出の性能が良いのかを理解するために、BatMeth2が呼び出したメチル化部位をResult A、BatMeth2-noIndelとBismarkが呼び出したメチル化部位をResult Bと定義してみました。 その結果、mclAには23,853個のDNAメチル化部位が含まれ、BatMeth2が呼び出したインデルリードをインデル検出でアラインメントした結果、23,853個のうち15,048個（63％）がカバーされていた（追加ファイル1：図S3Cを参照）。 また、Result AとResult Bのインデル率は、それぞれ5%と0%に過ぎないことがわかった。 このことから、正確なインデル検出はDNAメチル化判定を向上させることができると結論づけられました。
DNAメチル化データの可視化
BatMeth2には、メチル化データを可視化するためのツールが用意されています。 BatMeth2の可視化機能を説明するために、（1）ヒトH9細胞株から117Gbasesのシングルエンドリード、（2）ヒトIMR90細胞株から105.2Gbasesのシングルエンドリード、（3）野生型イネから12.6Gbasesのペアエンドリードをダウンロードした。 まず、BatMeth2は染色体レベルでのシトシンメチル化密度を可視化することができます。 図4aの点は、100kbのスライディングウィンドウを50kbのステップで表しています。 また、個々のCpGあるいは非CpG部位におけるMLをゲノムブラウザで閲覧できるように、bedおよびbigWig形式のファイルも提供している（図4b）。 遺伝子やTEの密度と比較したところ、MLはTE密度と相関があり、遺伝子密度とは逆相関であることが観察された（図4c）。 このような傾向は、以前にもイネで観察されている。
染色体スケールにおけるメチル化レベルの可視化 a ヒト第10染色体におけるメチルシトシン密度を示す。 100Kbを50Kb単位としたスライディングウィンドウで、メチル化レベルを表す。 b ヒト第10染色体におけるH9細胞株とIMR90細胞株のDNAメチル化レベルおよびメチル化差領域（DMR）の分布の例。 c イネゲノム全体における遺伝子、トランスポゾン要素（TE）の密度およびDNAメチル化レベル。 パネルAはBatmeth2から生成された結果です。 パネルBは、Batmeth2
のBEDファイルを用いてUCSCブラウザで可視化した結果である。 より正確には、BatMeth2は、遺伝子の2kb上流、転写開始部位（TSS）、遺伝子本体、転写終了部位（TES）、遺伝子本体の2kb下流のMLを可視化することが可能である。 上流、本体、下流を比較すると、Fig.5aは、遺伝子本体のDNA MLがプロモーター領域より高いことを示している。 また、5つの領域を比較すると、明らかにTSS領域に谷があることがわかる（Fig.5b）。 BatMeth2は、イントロン、エクソン、遺伝子間領域、TE周辺のMLプロファイルを計算することもできます（追加ファイル1． Figure S4）。 さらに、BatMeth2は、異なるサンプルの全体的な遺伝子MLを便利に比較するために、遺伝子領域ごとに複数の遺伝子のヒートマップを提供できる（図5c）
異なるコンテキストでのDNAメチル化の可視化。 a 遺伝子上流、遺伝子本体、遺伝子本体下流の2Kb領域におけるDNAメチル化レベル b 遺伝子間のDNAメチル化レベルの凝集プロファイル c 遺伝子上流、遺伝子本体、遺伝子本体下流の2Kb領域における全遺伝子のヒートマップ
3番目に、DNAメチル化の分布を視覚化できるBatMeth2が挙げられた。 Additional file 1: 図S5Aは、H9とIMR90の細胞株におけるDNAメチル化分布を示したものである。 図中、DNA MLは、メチル化（M：> 80％）、部分メチル化とメチル化の中間（Mh：60-80％）、部分メチル化（H：40-60％）、非メチル化と部分メチル化の中間（hU：20-40％）、非メチル化（U： < 20％）の5つに区分けされる。 Additional file 1に示すように 図S5Aに示すように、H9細胞株ではIMR90細胞株よりもMカテゴリでMLが高く、特にCpGコンテキストで高かった。 CH配列コンテキストでは、CpGメチル化が優勢であるが、CpA部位にはかなりの割合のメチル化シトシンが存在し、MLは特にH9細胞株で40％以下である（追加ファイル1：図S5B）
第四に、BatMeth2は遺伝子発現量と遺伝子プロモーターDNA ML間の相関を解析することができます。 H9とIMR90の細胞株を用いて、この機能を説明した。 H9やIMR90の遺伝子の発現量をカテゴリー別に分けてみた。 Additional file 1.に示すように 図S5Cに示すように、高発現遺伝子は、そのプロモーター領域で低いMLを示した。 さらに、遺伝子のプロモーターのMLを5つのカテゴリーに分類した。 その結果、Additional file 1: 図S5Dより、プロモーターのML値が高い遺伝子は、発現レベルが低いことがわかった。 哺乳類遺伝子の発現とプロモーターのDNAメチル化には負の相関があることが知られています。 この解析結果は、BatMeth2の精度の高さを示しています。
Finding Differentially methylated cytosines and regions (DMCs/DMRs)
メチル化データ解析において、DMCやDMRの同定は大きな目標の1つであります。 研究者は、単一のシトシン部位を表現型に相関させることに興味を持つこともありますが、DMRは非常に重要な特徴です。
初期のBS-Seq研究では、複製を収集せずに細胞をプロファイルしました。 このようなデータセットでは、フィッシャーの正確検定を使用して、メチル化シトシンの差（DMC）を識別した。 複製があるBS-Seqデータセットでは、DMCを呼び出すための最も自然な統計モデルは、ベータ二項分布です。 我々は、methylKit（生物学的複製を必要とする差分解析プログラム）やMethy-Pipe（生物学的複製を必要としない差分解析プログラム）など、多くのソフトウェアでDNAメチル化データの差分解析ができることを知っています。 しかし、マッピングとディファレンシャルメチル化解析の両方を含む包括的なパッケージは存在しない。 そこで、マッピングとディファレンシャル解析を統合したパッケージを開発した。 複製を持たないバイサルファイトデータからのDMRの同定を容易にするため、仮説検定を行うためにフィッシャーの正確検定を統合した。 サンプルが2つ以上の複製を持つ場合、ベータ二項分布を用いて差分メチル化解析を行う。 また、DMRのリストとしてbedファイルまたはbigWigファイルを提供しています。 図6aは、BatMeth2によって検出されたIMR90細胞株とH9細胞株におけるDMCと領域の数（p値< 0.05, meth.diff > = 0.6）を示したもので、DMRはゲノムブラウザ（図4b）上で可視化することができる。 BatMeth2は、遺伝子、CDS、イントロン、遺伝子間、UTR、TE、LTR、LINE、SINE領域など、ある領域でCpGやDMCが濃縮されているかどうかを可視化することが可能である。 図6bは、異なるゲノム領域におけるDMCの比率を可視化したものである。 遺伝子間領域を除けば、どの領域でもDMCの濃縮は見られなかった。
差分メチル化解析。 a H9細胞株とIMR90細胞株の差メチル化領域（DMR）、差メチル化遺伝子（DMGs）、差メチル化プロモーター（DMPs）の解析結果。 b 異なるゲノム特性および繰り返し要素に対するdifferentially-methylated Cytosine (DMC)のアノテーション。 c H9やIMR90に特異的なインデルを含むDMP（オレンジ色）は、全DMP（DNA Methylation differential Promoters）の中でかなりの割合を占める
A considerable proportion of differentially methylated promoters (DMPs) contains indels
組織形成や疾病においてインデルとDNAメチル化は、重要な役割をしていると知られています。 ここでは、differentially methylated promoters (DMPs)とインデルの関係について検討した。 IMR90とH9細胞株のBS-Seqリードを使用してこの研究を行った。 まず、BatMeth2を用いてBS-Seqリードをアライメントし、BisSNPとGATKツールを用いてインデルをコールした。 その後、H9またはIMR90のみに生じるインデルを細胞株特異的なインデルと定義した。
次に、BatMeth2によりH9とIMR90間で1384のDMPを検出した（p値< 0.05, meth.diff > = 0.6 ）。 上記の全DMPのうち、図6cに示すように、合計236個（17％）のDMPにインデルが含まれている。 つまり、かなりの割合のDMPがインデルを含んでいることがわかる。 したがって、これらのインデルの近傍のBS-Seqリードを正確にアライメントすることは、DNAメチル化の研究・探索に非常に重要である
。