macでインフォマティクス

macでインフォマティクス

NGS関連のインフォマティクス情報についてまとめています。

トランスポゾン検出ツール1 MELT

 

MELTは、iiluminaのペアエンドデータを使いリファレンスに存在しないmobile elementを検出するツール。以前1000 genomeで使われていたが、その後バージョンアップにより様々なゲノムに対応するようになった。SGEの分散コンピュータ環境から、SGEを使わない環境でも動く高いスケーラビリティを備えている。依存するソフトウエアも最小限になっており、エラーが起こりにくい設計を特徴とするらしい。

 

 

原理としては、一般的な構造変化の検出法と同じものを用いているが、途中からトランスポゾン検出に特化したパイプラインになっている。流れは大きく4つに分けることができる。1、ペアリードのマッピング(--> <--)が異常なペア(discordant read pair)の探索。2、おかしペアリードをユーザーが与えたトランスポゾン配列にアライメント。3、クロモソーム上を"waiking"して、トランスポゾンの挿入位置を正確に求める。4、情報を集計し、トランスポゾンごとに出力。この2のアライメントでbowtie2を使うため、MELTをランするにはbowtie2が前もってインストールされている必要がある。

 

ダウンロード

公式サイトのdownloadタブからダウンロードできる。

MELT - Home

解凍して中に入る。MELT.jarが実行ファイルである。

 

依存するもの

  • bowtie2
  • java -version1.7

 

2017年6月現在、javaはversion1.8が最新になっている。ランするにはjava 1.7が必要になってくるので注意してください(java7リンク)。=> その後1.8でも動作はしました。

 

ラン

入力はbwa memかalnで作成されたアライメントファイルである。  結果は標準的なvcf4.1のフォーマットで出力される。VCFという変異コールの標準的なフォーマットを採用しているので、MELTで解析後、VCFを扱えるツールに持ち込むことも容易になっている。ランは

java -jar MELT.jar

で行うが、いくつかオプションが必要である。

マニュアルに記載されているヒトゲノムのデータでテストランしてみよう。ラン前にhg19のbamとリファンレスゲノムをダウンロードする。

wget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/NA12878/alignment/NA12878.mapped.ILLUMINA.bwa.CEU.low_coverage.20121211.bam

bamのindexもダウンロードする。

wget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/NA12878/alignment/NA12878.mapped.ILLUMINA.bwa.CEU.low_

リファレンスゲノムをダウンロードする。

wget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz.fai

解凍する。

gunzip hs37d5.fa.gz

リファレンスのindexをダウンロードする。

wget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz.fai

リネームする。

mv hs37d5.fa.gz.fai hs37d5.fa.fai

トランスポゾンの配列を準備する。テストデータにはALU、SVA、LINE1の配列が用意されている。その名前をlsコマンドで出力。

ls /full/path/to/MELTvX.X/me_refs/1KGP_Hg19/*.zip > mei_list.txt

 exampleのトランスポゾン配列はzipで圧縮されている。

f:id:kazumaxneo:20170702132443j:plain

fastaで配列は記載され、bowtie2 samto olsでindexをつけられている。

f:id:kazumaxneo:20170702132557j:plain

 

 

準備ができたらランする。

java –jar MELT.jar Single -a –c 8 –h hs37d5.fa –l NA12878.mapped.ILLUMINA.bwa.CEU.low_coverage.20121211.bam –n add_bed_files/1KGP_Hg19/hg19.genes.bed –t mei_list.txt –w working_dir/

--single シングルゲノムの解析

perlで構築されているため、実行速度は遅めである。hg19の予測では、オーバーナイト走っていた。

指定したフォルダに結果は出力される。今回はALUS、SVA、LINE1の配列を指定したので、3つの配列の挿入部位がVCF4.1形式でコールされた。中身を確認する。

f:id:kazumaxneo:20170702113633j:plain

後半のフィールドが長くて見えづらい。前半のみ表示

SVA_final.vcf

f:id:kazumaxneo:20170702113751j:plain

LINE1_final.vcf

f:id:kazumaxneo:20170702114028j:plain

ALUは442コールあった。多すぎるので省略する。

 

 

トランスポソンタギングしたような株ではどうなるのかもテストしてみる。

検証中

 

 

注意点

トランスポゾン配列はfastaで指定する。名前はなんでもよいが、空白は許されていない。塩基配列は小文字、大文字とも認識可能で、60文字ごとに改行されていても1行に全部書かれていても動作する。ただしNがあってはいけない。MELTはATGCのみ認識する。また、トランスポゾン配列は1つずつファイルにする必要がある。マルチファスタ形式は許されていない。

  • BEDファイル

repeatのBEDファイルが必要である。そのため、UCSCのcustom track機能で

frepeat masker(11) trackをつけてダウンロードする方法が推奨されている。ヒトゲノムならば、UCSCのtable browser(Table Browser)から、GROUPS => Repeats、track => RepeatMaskerにチェックをつけてダウンロードする。

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

引用

An integrated map of structural variation in 2,504 human genomes.

Nature. 2015 Oct 1;526(7571):75-81. doi: 10.1038/nature15394.

Sudmant PH1, Rausch T2, Gardner EJ3, Handsaker RE4,5, Abyzov A6, Huddleston J1,7, Zhang Y8,9, Ye K10,11, Jun G12,13, Hsi-Yang Fritz M2, Konkel MK14, Malhotra A15, Stütz AM2, Shi X16, Paolo Casale F17, Chen J8,18, Hormozdiari F1, Dayama G19, Chen K20, Malig M1, Chaisson MJ1, Walter K21, Meiers S2, Kashin S4,5, Garrison E22, Auton A23, Lam HY24, Jasmine Mu X8,25, Alkan C26, Antaki D27, Bae T6, Cerveira E15, Chines P28, Chong Z20, Clarke L17, Dal E26, Ding L10,11,29,30, Emery S31, Fan X20, Gujral M27, Kahveci F26, Kidd JM12,31, Kong Y23, Lameijer EW32, McCarthy S21, Flicek P17, Gibbs RA33, Marth G22, Mason CE34,35, Menelaou A36,37, Muzny DM38, Nelson BJ1, Noor A27, Parrish NF39, Pendleton M38, Quitadamo A16, Raeder B2, Schadt EE38, Romanovitch M15, Schlattl A2, Sebra R38, Shabalin AA40, Untergasser A2,41, Walker JA14, Wang M33, Yu F33, Zhang C15, Zhang J8,9, Zheng-Bradley X17, Zhou W20, Zichner T2, Sebat J27, Batzer MA14, McCarroll SA4,5; 1000 Genomes Project Consortium, Mills RE19,31, Gerstein MB8,9,42, Bashir A38, Stegle O17, Devine SE3, Lee C15,43, Eichler EE1,7, Korbel JO2,17.

https://www.ncbi.nlm.nih.gov/pubmed/26432246?dopt=Abstract