macでインフォマティクス

macでインフォマティクス

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

リファレンスフリーで家族内変異や病変組織の変異を調べ、数十以下まで候補を絞り込む DIAMUND

追記 4/16

エラーが大量に出たので内容を修正しました。

 

   遺伝性疾患と癌の両方を含む、疾患の原因である突然変異を発見するためのゲノムシーケンシングの使用は、近年爆発的に増加している。全ゲノムシーケンスおよび全exomeシーケンスは、疾患表現型の原因となる可能性のある遺伝的変化を検出するために、何千人もの個体に採用されている。多くの成功は、ミラー症候群(Ng et al、2010)(わずか4人の脾臓の配列決定後に発見された)、バーター症候群[Choi et al、2009]の突然変異を同定するためのexomeシーケンシングを使用した2009年の報告に続き、フリーマン - シェルドン症候群[Ng et al、2009]に続いている。

   現在、高品質でほぼ完全なGRC37 / hg19ヒトリファレンスゲノム[The International Human Genome Sequencing Consortium、2001]は、これらの研究のすべてにおいて重要な役割を果たしている。分析のための標準プロトコールは、発端者および家族からのすべての配列をGRC37にアライメントさせ、次いで、一塩基多型(SNP)、小さな挿入および欠失(indels)、およびコピー数バリアントを検出する。しかし、ヒト個体群の自然変動のために、無作為の個体とGRC37とを比較すると、エクソームのみでも5万〜10万個の変異が生じる。ゲノム全体のシーケンスでは、数百万もの変異をもたらす。例えば、Roach et al(2010)は、ファミリーカルテットとヒトリファレンスゲノムを比較し450万個のSNPを見つけている。この膨大な候補の中で、しばしば単一の突然変異であるメンデリア病を引き起こす変異を見出すことは、乾草の中で針を探す (needle in a haystack) に類似したものになってしまう。最近の研究 (論文執筆時点) では、Strelka [Saunders et al、2012]およびFamSeq [Peng et al、2013]のような候補の膨大なセットを減らす統計モデルの開発に焦点が当てられてきたが、本質的な問題は残っている。リファレンスゲノムへのマッピングによって解析を始めると、非常に多数のリファレンスとの違いを検出することになる。

   Exome分析アルゴリズムでは、一連のフィルタを使用して50,000種類以上のバリアントのセットを管理可能なセットに絞り込み、個別に評価し、フォローアップ実験で検証することができる。しかしこれらのフィルタは、時には対象となる真のバリアントを除外する。例えば、フィルタリングの基準では、dbSNPのような大きなSNPデータベース、または1000Genomes Projectのような大きなプロジェクトで頻度> 1%で観察される変異すべてを除外する。これはフィルタとして非常に効果的であり、共通の母集団バリアントを表す何千もの無害なSNPを除去できる。しかし、このステップは、これらのデータベースのSNPのどれもが病気を引き起こすものではないことを暗黙のうちに仮定している。不完全な浸透、被験者の認識不能な疾患、およびその他の要因のために、健常者のみからSNPが収集されたとしても、この仮定は誤っている可能性がある。 解析パイプラインは、通常、多数の変異を除去するが、機能的に重要であるかもしれないスプライス部位および非コード変異も排除する。これらのフィルタや他のフィルタは意図せずに真の病因を排除する恐れがあり、さらにフィルタリング処理後でもexomeシーケンスで見つかる突然変異の数はまだ圧倒的に多いのが普通である。この問題は全ゲノムシークエンシングではさらに悪化し、無関係な変異の数は50倍になる。

   この論文では、exomeおよび全ゲノム解析に異なるアプローチをとるDiamund (direct alignment for mutation discovery) という新しい方法を紹介する。この方法では、突然変異候補は大幅に少なくなる。全てのサンプルをリファレンスゲノムにマッピングさせるのではなく、リードを互いに直接アライメントさせる。この方法は、(1)病変組織を同一個体の正常組織と比較する自己比較、および(2)被験体からのDNA配列間の差異がリファレンスとの比較よりずっと少ない家族内での解析で特に上手く機能する。 この方法では、germlineのリード (一般的にexome全体で1億以上の数)をGRC37リファレンスゲノムに合わせる必要はなく、複雑なゲノムアセンブリやこれらの大規模データセットからの統計的フィルタリングも使わない。以下で詳細に説明するように、著者らはどのサンプルでもユニークなシーケンスを素早く見つけることができる、より効率的なアルゴリズムを使用する。著者らはDiamundを2種類のexomes解析でテストした。まず、罹患した個体の病原組織に由来する初代培養線維芽細胞からのDNAを、同一個体の非病原性の初代培養線維芽細胞のDNAと比較したデータを検討した。腫瘍細胞または他の体細胞モザイクの遺伝的異常の分析のために、この直接的な比較は、最初に全てのリードをリファレンスゲノムと比較する分析よりも小さな変異数になる。第二に、family triosのデータを使った。このケースでは、子供のde novo突然変異が病気を引き起こしている疑いがあった。標準的なアルゴリズムでは、3つ全ての個体をリファレンスゲノムと比較するため、非常に大きな変異リストが出来上がるが、その多くは子供と親によって共有されるはずである。子供のDNAを両親と直接比較することにより、感度を失うことなく、またプロセスにノイズを加える家族特有の共通変異を検出することなく、すべてのde novo突然変異を迅速に同定することができた。真の変異を除いてしまう恐れのあるフィルターに頼る方法を取り除くことで、これらのテストでは、真のde novo変異リストは非常に少なかった。

   De novo突然変異は、メンデル症の高い割合を占める可能性がある。 Yang et alは250人の発端者(probands)およびその家族のexomeシークエンスについて調べ、常染色体優性患者33人およびX連鎖病の9人を同定した。これらのうち、常染色体優性の83%およびX連鎖変異の40%が新たに発生したものだった。ファミリー内のサンプル間、または影響を受けた組織と影響を受けていない組織との間の直接的な比較は、偽陽性を減らし、リファレンスゲノムから完全に欠けている領域の突然変異の検出も可能にする。いくつかのヒト集団は、多くのメガベースにまたがる大きなゲノム領域を共有していることが既に示されている[Li et al、2010]。これはヒトリファレンスゲノムから完全に欠けている。これらには、新規の断片的重複[Schuster et al、2010]ならびに全く新しい配列が含まれる。関心のある突然変異がこれらの領域の1つに落ちた場合、従来の方法はそれを見逃す。対照的に、著者らの直接比較アルゴリズムは、これらの領域も除かないため、それらの領域内での突然変異を見つけることができる。 重要な注意点として、Diamundは変種検出のよりgeneralな問題を解決することを意図していない。この方法論は、サンプルの1つまたはサブセットに存在する突然変異をより効果的に同定することができうる、非常に密接に関連するサンプルを利用する前提で設計されている。

 

Method

 2つ以上のゲノムを直接比較する1つの方法は、データセットをde novoでアセンブルし、全ゲノムアライメントが可能なMUMmerなどで比較することである。しかし全ゲノムアセンブリは計算コストが高く、かつ誤ったアセンブリを生成する可能性もあり、アライメントによるアプローチよりもさらに大きな問題を引き起こす可能性がある。代わりにDiamundは、ある固定サイズのkについて、すべてのリードで長さkの全シーケンスを数え、次にこれらのk-mersを互いに比較する直接的なアプローチを使用する。アルゴリズムの10の主要なステップを概説する。

 

1、jellyfishを使ったk-mer抽出と頻度カウント。シーケンスしたリードが全て正しければ、(偶発的なマッチがほぼ起こらない)大きめのサイズでk-merをカウントすると、異なるk-merの数はゲノムのサイズ(ここではexomeのサイズである50-60Mb) と等しくなる (ゲノムにリピート領域がない時)。しかし実際は、シーケンスエラーにより異なるk- merの数は爆発的に増える。これにstep2で対処する。

 ↓ 

2、シーケンスエラー由来k-merのフィルタリング。リードに1つのシーケンスエラーが発生すると、本来ゲノムに無いだろうk個の新しいk-mer配列が出来る。そのようなk-merの出現回数は大抵1だが、exome解析ではカバレッジは100近くあり、偶然2以上の出現頻度になるシーケンスエラー由来のk-merは、ゲノム全体では無視できない数になる。著者らのempericalな観測から、頻度10以下のk-merを除くことで、真の情報の損失は限りなくゼロにしつつ、シーケンスエラーはほぼ完全に排除できると判断した。頻度10以下のk-merを捨てる。

  ↓ 

3、ベクターとリファレンスヒトゲノムに出てくるk-merのフィルタリング。ヒトリファレンスゲノムのk-merを除く。これは、リファレンスゲノムに存在するk-mer配列は病気の発生に関与しないと考えられるためである。ベクターデータはUniVec database を使っている。

  ↓ 

4、Proband (子供) と父親、probandと母親それぞれの共通/非共通のk-merの探索。両親のシーケンスデータには見つからないk-merを抽出し、父との比較、母との比較、得られた二つの結果をミックスする。得られたk-merには、probandのexomeに固有のk-merが全て含まれているはずである。

  ↓ 

5、固有のk-merをリードにバックアライメントして、リードに差し戻す。ステップ4の終了時点で情報は1万倍以上濃縮されており、論文のデータでは、数千〜数万のk-merまで減っている(論文 表1)。このk-merを含むリードを抽出し、リードに戻す。MUMmerかmodifyしたkrakenを使いシーケンスリードとのアライメントを行うことで、該当リードを見つけ出す。デフォルトではより高速なkrakenを使う。

  ↓ 

6、再度ベクター配列のフィルタリング。probandsにまだ残っているベクターデータを除く。ステップ3と同じデータベースを使うが、精度の高いvecscreen programを使っている。初期インプットと比べリード数は極端に減っているので、この段階では、コンピュータリソース負荷が大きいプログラムの使用は容易になっている。

 

7、選抜したリードのアセンブリMinimoアセンブラを使い、リードをアセンブルする。ここまでの処理で残ったリードは、変異を含む部位に限定されているはずで、アセンブルしても伸びることはない(変異が集中しない限り)。100-bpのリードならば、アセンブルしても、理屈上、最大199-bpになる。実際はもっと短くなる。

  ↓ 

8、コンティグの選抜。シーケンスエラーにより、同じポジションでもリードにバリエーションが生じ、複数コンティグが生じることがある。bowtieを使ってコンティグをリファレンスゲノムにマッピングする。ゲノムの同じ部位に複数のコンティグがアライメントされるならば、一つのコンティグに統合する。100-bp以下の配列は、100-bpまでリファレンス情報を使い拡張する。(リファレンスゲノムはあくまでコンティグのクラスタリングのためにしか利用していない)。

  ↓ 

9、コンティグの選抜2。両親の両方または片方のexome キャプチャに失敗して、probandに固有の配列に見えている可能性がある。bowtieを使って両親のシーケンスリードを8のコンティグにマッピングして、両親のデータともマッピングされるか確かめる(この時、真のマッピング部位か検証するため、マッピングされたリードは、リファレンスゲノムにもマッピングし、ベストマッピングかどうか確認している)。片親データしかマッピングされなければ、実験のエラーの可能性が捨てきれないので排除する。

またstep8でリファレンスゲノムにアライメントされなかったcontigは、別に保存される(情報がないのでマニュアルで解析する必要がある)。

  ↓ 

10、バリアントのコールとアノテーションsamtoolsのmpileupを使い、リファレンスゲノムと比較して変異している部位を、アミノ酸変化の注釈付きで出力する。

f:id:kazumaxneo:20180410123520j:plain

論文より転載。

 

 

公式サイト

http://ccb.jhu.edu/software/diamund/ 

インストール 

cent OS6に導入した。

 依存

  • Kraken a system for assigning taxonomic labels to short DNA sequences;
  • Samtools package for managing large short read alignment files;
  • Jellyfish a tool for fast, memory-efficient counting of k-mers in DNA;
  • Bowtie2 an ultrafast and memory-efficient tool for aligning sequencing reads;
  • Minimo a de novo assembler based on the AMOS infrastructure;
  • Mummer a system for rapidly aligning entire genomes.
  • vecscreen
  1. Minioは上記リンクからダウンロードして、"./configure"、"sudo make install"する(binにMinimoのバイナリができる)perlのライブラリも必要になるので、cpanm XML::Parserで導入する(DBIStatistics::Descriptiveも同様)。
  2. krakenとjellyfishは以前の説明を参照(リンク)。jellyfishはversion1をインストールすることに注意する。version1ダウンロード方法はkrakenの解説を参照。
  3. 他のツールはbrewで導入できる。
  4. <追記> SAMToolsのコードが旧式のものになっている。旧式のSAMtoolsバイナリをダウンロードして(リンク)、そのSAMtoolsを指定する。 mpileupを使うので、version 0.1.13あたりを使う。perlのコードを現在のSAMtoolsの形式に修正してもよいが、いくつかのサブルーチン内をいじるので手間になる。ソースをダウンロードしてmakeすればsamtoolsができるので、ここにdiamond.cfgのsamtoolsパスを通す。現行のsamtoolsと使い方が少しだが致命的に異なるコマンドがあるので、旧式を/usr/local/bin/samtoolsなどと置き換えてはいけない(sortのコマンドとかmpileupのオプションの違いとか)。
  5. <追記> vecscrrenは公式には記載がないが必要になる。ダウンロードできるURLが見つからなかったので、Githubのvecscreen_plus_taxonomyレポジトリのscritsにあるvecscreenのバイナリをダウンロードした(linux版)(ただしこのvecscreenはバージョンが違うので、Diamondのコードを修正する必要がある)。

    vecscreen_plus_taxonomy/scripts at master · aaschaffer/vecscreen_plus_taxonomy · GitHub

 

本体は公式からダウンロードして解凍する。

tar -xvf diamund.tar
cd Diamund/
perl diamund.pl 

$ perl diamund.pl 

Usage:

 diamund.pl {<patient>.fa|<patient>.fq[.gz|.gz2][,<patient_mates>.fq[.gz|.gz2]]} {<parent1>.fa|<parent1>.fq[.gz|.gz2][,<parent1_mates>.fq[.gz|.gz2]]}  [{<parent2>.fa|<parent2>.fq[.gz|.gz2][,<parent2_mates>.fq[.gz|.gz2]]}] {-e <exon_file>.fa|-x <exon_kmer_file>} [options]

Options:

 -e <exon_file>         fasta file with the sequence of all exons

 -x <exon_kmer_file>    sorted exome kmer file (either the -e or -x options should be present in the input)

 -k|kmer <kmer_size>    kmer size for jellyfish program (default: 27)

 -r|recessive           run program in recessive mode (default: dominant)

 -t <no_of_threads>     number of threads for the jellyfish program

 -v <kmer_vector_file>  file with vector kmers (if available)

 -z <vector_file>       fasta file with vectors

 -f <filter_threshold>  only keep kmers in patient that appear more than the <filter_threshold> times

 -n <read_no>           only keep assemblies that have at least <read_no> assembled reads 

 -d <vecscreen_data>    screen for vector contamination with vescreen using the data from this directory

 -g <[0|1](,[0|1])*>    a string of 0 and 1 separated by commas, and signifying the sex associated with each data sample respectively; 1 is male, and 0 is female  

 -l <read_len>          read length (default:100)

 -o <output_file>       output file with variants called (default is STDOUT)

 -s <sequence_depth>    observed depth of coverage in non-patients in order to call variant (default:10)

 -y <kmer_no>           assume kmers occuring this many times in normals are errors (default: 0)

 -c                     don\'t clean-up the files                  

 -m                     use mummer instead of kraken to align kmers to reads (default: no)

Usage examples:

 diamund.pl patient.fa father.fa mother.fa -e hg19_exome.fa -g 1,0,1

or

 diamund.pl patient.1.fq.gz,patient.2.fq.gz normal.1.fg.gz,normal.2.fg.gz -x 27mers-from-exome-sort

Error: at least two input fasta file need to be specified

 

 公式からhuman hg19のFASTAとbowtie indexをダウンロードする。

wget ftp://ftp.ccb.jhu.edu/pub/software/diamund/data/hg19.tar.gz
mkdir human_ref_index
tar -zxvf hg19.tar.gz -C human_ref_index/ #human_ref_index/に解凍

Vectorとexomeの計算済みk-merデータベースも公式に準備されている。ダウンロードする。

wget ftp://ftp.ccb.jhu.edu/pub/software/diamund/data/27mers_exome_vect.tar.gz
mkdir k-mer_databae
tar -zxvf 27mers_exome_vect.tar.gz -C k-mer_databae/

ベクターVeccscreenのデータベース配列へのリンクも公式に載っているが、NCBI ftpサーバーのURLが変更されている(The UniVec Database)。以下のURLでダウンロードできる。

wget ftp://ftp.ncbi.nlm.nih.gov/pub/UniVec/UniVec

UniVecはFASTAフォーマットのデータ。

$ head UniVec

>gnl|uv|X66730.1:1-2687-49 B.bronchiseptica plasmid pBBR1 genes for mobilization and replication

CTCGGGCCGTCTCTTGGGCTTGATCGGCCTTCTTGCGCATCTCACGCGCTCCTGCGGCGGCCTGTAGGGC

AGGCTCATACCCCTGCCGAACCGCTTTTGTCAGCCGGTCGGCCACGGCTTCCGGCGTCTCAACGCGCTTT

GAGATTCCCAGCTTTTCGGCCAATCCCTGCGGTGCATAGGCGCGTGGCTCGACCGCTTGCGGGCTGATGG

TGACGTGGCCCACTGGTGGCCGCTCCAGGGCCTCGTAGAACGCCTGAATGCGCGTGTGACGTGCCTTGCT

GCCCTCGATGCCCCGTTGCAGCCCTAGATCGGCCACAGCGGCCGCAAACGTGGTCTGGTCGCGGGTCATC

TGCGCTTTGTTGCCGATGAACTCCTTGGCCGACAGCCTGCCGTCCTGCGTCAGCGGCACCACGAACGCGG

TCATGTGCGGGCTGGTTTCGTCACGGTGGATGCTGGCCGTCACGATGCGATCCGCCCCGTACTTGTCCGC

CAGCCACTTGTGCGCCTTCTCGAAGAACGCCGCCTGCTGTTCTTGGCTGGCCGACTTCCACCATTCCGGG

blastのデータベースを作る。

makdir database
makeblastdb -dbtype nucl -in UniVec -out database/UniVec
ls database/

$ ls database/

UniVec.nhr  UniVec.nin  UniVec.nsq

 

 diamond.cfgを開き、各ツールのパスを記載しておく。さきほどダウンロードしたデータベースのパスと、コマンドのパスから(whichで調べる)、以下のように変更した。

$ cat diamund.cfg 

mummer  ~/.linuxbrew/bin/mummer

vecscreen_data  /home/uesaka/Diamund/database/UniVec

kraken  /home/disk2/uesaka/kraken

human_index /home/uesaka/Diamund/human_ref_index/hg19

human_fasta /home/uesaka/Diamund/hg19.fa

jellyfish /home/disk2/uesaka/kraken/jellyfish-1.1.11/bin/jellyfish

minimo  ~/amos-3.1.0/bin/Minimo

bowtie  ~/.linuxbrew/bin/bowtie

bowtie2  ~/.linuxbrew/bin/bowtie2

samtools /home/uesaka/Diamund/samtools-0.1.13/samtools

bcftools  /usr/local/bin/a5_miseq_linux_20160825/bin/bcftools

vecscreen /user/local/bin/vecscreen

 bowtieのほか、はbowtie2の行を追加する。 

 

vecscrren2を使ったので、オプションを変更する必要があった。結果が変わる可能性を覚悟で988行目の

my $out= `echo $seq | $vecscreen -d $vecscreen_data -f3`;

my $out= `echo $seq | $vecscreen -db $vecscreen_data`;

に変更した。

 

 

 

 

ラン

 fammily triosの解析。3人のシーケンスデータを入力する。

diamund.pl patient.fa father.fa mother.fa -e hg19_exome.fa -g 1,0,1 
  • -e <exon_file>     fasta file with the sequence of all exons 
  • -g <[0|1](,[0|1])*>   a string of 0 and 1 separated by commas, and signifying the sex associated with each data sample respectively; 1 is male, and 0 is female

正常組織と病原組織の解析。fastqはgz圧縮形式にも対応している。

diamund.pl patient.1.fq.gz,patient.2.fq.gz normal.1.fg.gz,normal.2.fg.gz -x 27mers-from-exome-sort
  •  -x <exon_kmer_file>   sorted exome kmer file (either the -e or -x options should be present in the input)

jellyfishはメモリが足りないとsegmentation errorを起こす。WESでも複数走らせる時は数百GBいるかもしれない。注意する。

 

 

テストデータとして、 Manuel Corpasさんが公開してくれている親子トリオのexomeデータを使ってみる。

https://cambridgeprecisionmedicine.com/2011/07/12/getting-my-genome-sequencing-done-part-i/

https://cambridgeprecisionmedicine.com/2013/01/21/corpas-family-exome-data-available-for-public-download/amp/

slideshareからダウンロードし解凍する。

unzip 106340.zip

dad_1.clean.fq.gz、dad_2.clean.fq.gz、mom_1.clean.fq.gz、mom_2.clean.fq.gz、daughter_1.clean.fq.gz、daughter_2.clean.fq.gzができる。

父、母、娘のデータを指定してランする。

perl diamund.pl dad_1.clean.fq.gz,dad_2.clean.fq.gz mom_1.clean.fq.gz,mom_2.clean.fq.gz daughter_1.clean.fq.gz,daughter_2.clean.fq.gz -s 10 -t 40 -x k-mer_databae/27mers-exome

 

krakenを使ってk-merからリードに戻すstepでエラーになる。ツールのバージョン違いの気が下が、手っ取り早く-mをつけてkrakenからmummerに切り替えて対応した。

bcfしか出力しなかったので、bcftoolsを使いvcfに変更した。

bcftools view 2-to-1.bcf > out.vcf

個人データの結果は載せないが、出力は確かにかなり限定された。

 

 

 

感想

思いっきりempericalなアルゴリズムwiki)に頼った手法なので、条件を満たせば良い結果を出してきます。数学的な拠り所が少ないのでインフォマティクスが専門の方には興味が持たれないツールかもしれませんが、wetの研究者にとっては良いツールではないでしょうか?

 

ただし、このツールはどうやら堅牢性があまり考慮されておらず(本体はperlスクリプト)、ツールのバージョンが少しでもずれるとエラーを起こし、そこでジョブが止まります。マニュアル不足で3rd partyのツールのインストールに手間取ったのも影響してますが、繰り返しテストしていたら1週間以上かかりました。ランニング時は > log.txt 2>&1とかをつけて、標準出力のエラーログを取りながら、どのstepでつまづいたか1つ1つ確認して解決していくのが良いでしょう。またテスト時はchr1の一部などにマッピングして、マッピングされたfastqだけ回収し、それで進めてください。速度が考慮されておらず、テストにもかなり時間がかかります。

 

現在開発中となってますが、RUFUSも発想は同様のツールのようです(ただしユニークk-merを抽出してそれをそのままアセンブリする)。

GitHub - jandrewrfarrell/RUFUS: RUFUS k-mer based genomic variant detection

 

 

引用

DIAMUND: Direct Comparison of Genomes to Detect Mutations

Steven L Salzberg, Mihaela Pertea, Jill A Fahrner, and Nara Sobreira

Hum Mutat. 2014 Mar; 35(3): 283–288.