macでインフォマティクス

macでインフォマティクス

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

高感度なバリアントコーラー LoFreq

 

 シーケンシング技術の最近の進歩により、細胞集団におけるheterogeneityおよびsub-populationsのより広範な研究、およびそれらの進化による「コンセンサス配列」からの移行が可能になった。このような‘population perspective’ な視点は、ウイルス性疑似種(Viral quasispecies、wiki)およびウィルスの宿主内バリエーション(論文より ref.1,2)、バクテリアのsub-populations(ref.3-5)、がんのサブクローン進化(ref.6-8)など適用可能な様々なアプリケーションを持つ。これらの研究における集団構造(およびレアな亜集団)の正確な特徴付けは、宿主応答または抗生物質曝露のダイナミクスにとって基本的である。いくつかの最近のガンシーケンシング研究はレアなsub-populationsおよびバリアントの機能的役割を、腫瘍のgrowth抗生物質耐性、転移(9,10)およびそれらを研究するための計算ツールの観点でさらに強調している。

 原理的には、大規模並列シーケンシングの高スループットにより、レアなsub-populationsのサンプリングが可能になる。しかし、シークエンシングエラーはpopulationsにおける真の変異の決定を複雑にする。シークエンシングエラー率は、テクノロジー、ラン、レーン、マルチプレックス、ゲノムロケーション、および置換タイプ(ref.11〜13)により著しく異なることが知られている。これらの問題を解決するためのアプローチが研究されてきたが、バリアントコール法の大部分は、関心のある個別の頻度(すなわち0,0.5および1; 低カバレッジのヒトリシーケンシングデータおよび二倍体コール(ref.14-16)関連の一連の方法は、プールされたシーケンシングデータ(17-20)において二倍体genotypesをコールするために調整されたものであり、一般には適用されない)。

 RNAウイルスの研究はこの規則(ref.21-24)の例外を提供している。 RNAウイルスは、RNA依存性DNAポリメラーゼのプルーフリーディング能力が不十分または欠損しているために高い突然変異率および高い複製速度を有し、結果として生じる宿主内変異は抗生物質治療戦略(ref.25)、生ワクチンの遺伝子モニタリングにも影響がある(ref.26)。これらの研究で使用された手法は ad hoc なトリミング、フィルタリング、閾値を使用してバリアントコールに頼っており、そのため感度や広範な適用性が制限されている(サンプルまたはシーケンスごとの手動調整が必要)。 Breseq(ref.27,28)やSNVer(ref.29)などのモデルベースの手法は、潜在的により敏感で一般的だが、単純な二項モデルに依存しており、エラー率の偏りに対応するように調整されていない。シーケンシングエラーの影響を減らすフェージングに依存し、454シークエンシングに合わせたより洗練されたアプローチが、ウイルスデータセットに適用されている(ref.30)。しかし、この方法は他の技術に直接適用することはできず、大きなゲノムやシーケンシングデータセットでは実行できない。

 シーケンシングを用いて細胞のゲノム状態をモニタリングする新しい臨床応用では、population内のレアなバリアントを検出し、検出限界の末端までそれを行う能力は、未達成の重要な能力である。一方で、バリアントコールの感受性の増加は、レアだが重要なsub-populations(e.g. cancer stem cell mutations)をモニターすることを可能にする。他方、薬剤耐性sub-populationsの早期検出には感度が不可欠である(e.g. with antiretroviral drugs for HIV)。そのような設定では、ad hocなアプローチは適合性およびロバストさを欠いており、バリアント検出感度において人工的なキャップを被る可能性がある。シミュレーションエラーの正確なモデリングは、感度限界を押し上げるために不可欠であり、取り組むべき必要性がある。

 この研究では、シークエンシングエラー率の偏りの公式モデルに基づいて、高カバレッジシークエンシングデータセットからシングルヌクレオチドバリアント(SNV)をコールための、高感度で堅牢なアプローチを提示する。近似に頼ることなく、このモデルでは厳密な統計的検定を効率的に行うことができ、大きなゲノムや高カバレッジのデータセットを正確に解析することができる。結果のメソッドLoFreqは、シーケンス実行およびポジション特異的なシーケンスバイアスに自動適応し、データセットの平均シーケンシングエラーレートより低い頻度のSNVを呼び出すことができる。 LoFreqの堅牢性、感度および特異性は、いくつかのシミュレートされたデータセットと実際のデータセット(ウイルス、バクテリアおよびヒト)および2つの実験プラットフォーム(FluidigmおよびSequenom)を用いて検証された。レアなsomatic SNV(胃がんのエキソームシーケンシングデータセット)コールにLoFreqを適用し、臨床研究(ヌクレオシドアナログ薬Balapiravir)で治療前後のデングウイルス quasi species を研究した結果は、さらに堅牢性と多様性を強調している。

 

 LoFreq(バージョン2)は、次世代シーケンシングデータからSNVとindelを推論するための、高速で機密性の高いバリアントコーラー。これは、通常他の方法によって無視される、ベースコールクオリティおよびシーケンスに固有の他のエラー源(例えばindelアラインメントの不確実さ)を十分に利用する。LoFreq *は、シーケンシングマシンまたはシーケンシング技術に依存しない閾値を使うため、ほぼすべてのタイプのアライメントしたデータ(例えばIllumina、IonTorrentまたはPacbio)で使用できる。またカバレッジおよびシーケンスクオリティ変化に自動的に適応するので、様々なデータセットにも適用できる(例えばウイルス、バクテリア、メタゲノムまたは体細胞シーケンシングデータなど)。LoFreqは非常に敏感で、平均ベースコールクオリティ(すなわち、シーケンシングエラー率)を下回るバリエーションを予測できる。各バリアントコールには、false postiveの確率をコントロールできるp値が割り当てられる。近似やヒューリスティックは使用しない。(擬似)並列実装も行われている。カバレッジに自動適応するので、カバレッジの高いデータや大きなゲノムに利用できる。単一プロセッサーでは、600XカバレッジのE.coliゲノムのSNVコールでおよそ1時間、100Xカバレッジのヒトexomeデータセットではおよそ1時間かかる。

 

LoFreq · Fast and sensitive variant calling from next-gen sequencing data

f:id:kazumaxneo:20180817231611p:plain

マニュアル

Commands · LoFreq

LoFreqに関するツイート。


インストール

mac os10.13でテストした。

依存(ソースからコンパイルする場合に必要)

  • a C compiler (e.g. gcc or clang)
  • a Python 2.7 or Python 3 interpreter
  • zlib developer files
  • a compiled version of samtools 1.1
  • a compiled version of htslib 1.1; use the one that comes bundled with samtools!)

Github

SourceForgeにバイナリが用意されている。またcondaでもインストールできる。

https://sourceforge.net/projects/lofreq/files/

ダウンロードしてパスを通す。

$ lofreq

 

       |             ____|                 

       |       _ \   |     __|  _ \   _` | 

       |      (   |  __|  |     __/  (   | 

      _____| \___/  _|   _|   \___| \__, | 

                                        _| 

 

Fast and sensitive inference of SNVs and indels

 

Usage: lofreq <command> [options]

 

  Main Commands:

    call          : Call variants

    call-parallel : Call variants in parallel

    somatic       : Call somatic variants

 

  Preprocessing Commands

    viterbi       : Viterbi realignment

    indelqual     : Insert indel qualities

    alnqual       : Insert base and indel alignment qualities

 

  Other Commands:

    checkref      : Check that reference fasta and BAM file match

    filter        : Filter variants in VCF file

    uniq          : Test whether variants predicted in only one sample really are unique

    plpsummary    : Print pileup summary per position

    vcfset        : VCF set operations

    version       : Print version info

 

  Samtools Clones:

    faidx         : Create index for fasta file

    index         : Create index for BAM file

    idxstats      : Print stats for indexed BAM file

 

  Extra Tools (if installed):

    vcfplot       : Plot VCF statistics

    cluster       : Cluster variants in VCF file (supports legacy SNP format)

 

lofreq call -h

$ lofreq call -h

lofreq call: call variants from BAM file

 

Usage: lofreq call [options] in.bam

 

Options:

- Reference:

       -f | --ref FILE              Indexed reference fasta file (gzip supported) [null]

- Output:

       -o | --out FILE              Vcf output file [- = stdout]

- Regions:

       -r | --region STR            Limit calls to this region (chrom:start-end) [null]

       -l | --bed FILE              List of positions (chr pos) or regions (BED) [null]

- Base-call quality:

       -q | --min-bq INT            Skip any base with baseQ smaller than INT [6]

       -Q | --min-alt-bq INT        Skip alternate bases with baseQ smaller than INT [6]

       -R | --def-alt-bq INT        Overwrite baseQs of alternate bases (that passed bq filter) with this value (-1: use median ref-bq; 0: keep) [0]

       -j | --min-jq INT            Skip any base with joinedQ smaller than INT [0]

       -J | --min-alt-jq INT        Skip alternate bases with joinedQ smaller than INT [0]

       -K | --def-alt-jq INT        Overwrite joinedQs of alternate bases (that passed jq filter) with this value (-1: use median ref-bq; 0: keep) [0]

- Base-alignment (BAQ) and indel-aligment (IDAQ) qualities:

       -B | --no-baq                Disable use of base-alignment quality (BAQ)

       -A | --no-idaq               Don't use IDAQ values (NOT recommended under ANY circumstances other than debugging)

       -D | --del-baq               Delete pre-existing BAQ values, i.e. compute even if already present in BAM

       -e | --no-ext-baq            Use 'normal' BAQ (samtools default) instead of extended BAQ (both computed on the fly if not already present in lb tag)

- Mapping quality:

       -m | --min-mq INT            Skip reads with mapping quality smaller than INT [0]

       -M | --max-mq INT            Cap mapping quality at INT [255]

       -N | --no-mq                 Don't merge mapping quality in LoFreq's model

- Indels:

            --call-indels           Enable indel calls (note: preprocess your file to include indel alignment qualities!)

            --only-indels           Only call indels; no SNVs

- Source quality:

       -s | --src-qual              Enable computation of source quality

       -S | --ign-vcf FILE          Ignore variants in this vcf file for source quality computation. Multiple files can be given separated by commas

       -T | --def-nm-q INT          If >= 0, then replace non-match base qualities with this default value [-1]

- P-values:

       -a | --sig                   P-Value cutoff / significance level [0.010000]

       -b | --bonf                  Bonferroni factor. 'dynamic' (increase per actually performed test) or INT ['dynamic']

- Misc.:

       -C | --min-cov INT           Test only positions having at least this coverage [1]

                                    (note: without --no-default-filter default filters (incl. coverage) kick in after predictions are done)

       -d | --max-depth INT         Cap coverage at this depth [1000000]

            --illumina-1.3          Assume the quality is Illumina-1.3-1.7/ASCII+64 encoded

            --use-orphan            Count anomalous read pairs (i.e. where mate is not aligned properly)

            --plp-summary-only      No variant calling. Just output pileup summary per column

            --no-default-filter     Don't run default 'lofreq filter' automatically after calling variants

            --verbose               Be verbose

            --debug                 Enable debugging

 

 

ラン

 利用するbamは、false positiveを防ぐため、予めPCR duplicationを除き、indel周辺のリアライメントを行って、GATKのBQSR(ヒトでは既知バリアントデータベースがあるので利用可能)(公式)を行うことが推奨されている。

バリアントコール。

lofreq call -f ref.fa -o vars.vcf aln.bam
  • --call-indels    Enable indel calls (note: preprocess your file to include indel alignment qualities!)
  • -a    P-Value cutoff / significance level [0.010000]
  • -C   Test only positions having at least this coverage [1]
  • -d    Cap coverage at this depth [1000000]
  • --illumina-1.3    Assume the quality is Illumina-1.3-1.7/ASCII+64 encoded
  • --verbose    Be verbose
  • -S   Ignore variants in this vcf file for source quality computation. Multiple files can be given separated by commas

 

normalと比較してtumorのバリアントをコール。

lofreq somatic -n normal.bam -t tumor.bam -f hg19.fa \
--threads 8 -o output -d dbsnp.vcf.gz

 

高感度で高カバレッジなデータセットにも自動適応できることを利用して、バクテリアの研究室進化実験(シミュレーション)でのheterogenicなバリアントコールにも使われています。

https://www.ncbi.nlm.nih.gov/pubmed/28224054

追記

 引用
LoFreq: a sequence-quality aware, ultra-sensitive variant caller for uncovering cell-population heterogeneity from high-throughput sequencing datasets
Andreas Wilm, Pauline Poh Kim Aw, Denis Bertrand, Grace Hui Ting Yeo, Swee Hoe Ong, Chang Hua Wong, Chiea Chuen Khor, Rosemary Petric,  Martin Lloyd Hibberd, and Niranjan Nagarajan 

Nucleic Acids Res. 2012 Dec; 40(22): 11189–11201.

 

https://www.sciencedirect.com/science/article/pii/S2001037017300946

review article要約 ラージゲノムのシーケンシング解析

 

 

はじめに

 この記事はレビュー論文の要約です。チェックリスト、または思考を整頓するガイドとして使ってください。ただし、この要約で論文を読んだことにはなりません。時間が許す限り原著論文を読むことをお勧めします。review articleは各段落が一般論で構成されており、そのぶん読みやすくなっています。毎日読んでいれば綺麗な英語も身につくでしょう。

 

 

Abstract

 長い間、ゲノムシークエンシングプロジェクトはモデル生物に限定されており、実施するには大規模なコンソーシアムの協調的な努力が必要であった。ハイスループットシーケンシング技術の急速な進歩とバイオインフォマティクスツールの並行した発展は、この分野の民主化を促した。現在では選択した任意の生物のde novoドラフトゲノム構築が、生態と進化(eco‐evolutionary)および保全の共同体における個々の研究グループの手の届く範囲内にある。そのような目標のために費やされるかなりの努力とコストを考えると、重要な第一ステップは、手元の生物学的問題に対処するためにゲノム配列が必要かどうかを徹底的に検討することである。一旦この決定が下されれば、ゲノムプロジェクトはその生物およびゲノムドラフトの品質に関して慎重に計画せねばならない。ここでは、この分野の最先端技術を簡単に見直し、大規模かつ複雑なゲノムを中心にゲノムシーケンシング、アセンブリアノテーションに関するワークフローを段階的に紹介する。このチュートリアルでは、保全遺伝学の背景を持つ科学者を対象としているが、より一般的には、全ゲノムシーケンシングプロジェクトに携わる研究者にとって有用な実践的指針を提供する。

 

オーサーらは脊椎動物ゲノムのアセンブリを経験したバックグランドを持つ。この論文はNGSの包括的レビューではなく、スコープをlarge genomeのシーケンシング解析に限定している(RNA seqやバクテリアゲノムのアセンブリについては議論しない)。自然保護、種の保護の分野にもゲノムシーケンシング解析が広がっており、本論文は、その保護の観点で、モデル生物よりも非モデル生物のゲノム解析に重きを置いて記述されている。

 

本文要約

Basic considerations

 ゲノムアセンブリはリソース、時間、経験を必要とするチャレンジングな問題である。以下の点について慎重に検討しなければならない。

  1. そもそもゲノムシーケンシングが必要か?ゲノム配列は単なるリソースであり、conservation biologyなどの分野では役に立たないケースも多い。
  2. ゲノムアセンブリを行なっても、大きなファミリーを構成しているコピー遺伝子や、進化スピードの早い遺伝子は、最終的なドラフトアセンブリでもほとんど配列として表現されていない可能性がある。決定するには、別の手法を並行して使用したり、BACなどにクローニングして個別にシーケンシングする必要がある。
  3. ターゲットの領域が決まっているなら、そこだけPCR増幅して次世代シーケンス解析することで、低コスト且つ高いカバレッジでシーケンスすることができる。時間やコストも大きく削減できる(例 Wang et al., 2012)。Box 2. Before you startに考慮するべきこと一覧が挙げられている(box2 論文リンク

 

ゲノムをシーケンスするとはどうゆうことか?

 ドラフトゲノム配列には全クロモソームの全塩基配列が表現されているのが理想的であるが、実際には以下の理由からゲノムシーケンスのコンセプトに対して大きな矛盾、齟齬が存在している。

  1. ゲノム配列には個体ごとにバリエーションが存在する。よってある1つの種に対する1つの真のゲノム配列というものは存在しない。
  2. 2倍体(2n)生物のゲノムには、相同染色体間にヘテロ接合性のポジション、indel多型、コピー数変化、小規模なゲノムリアレンジメントが存在する。これも真のゲノム配列を定義することを困難にしている。
  3. 1つの個体の特定の細胞群だけシーケンシングしても、体細胞変異による後天的変異が一定の確率で起こり得るため、これもゲノムにバリエーションを引き起こす。
  4. 2倍体以上のポリプロイドの種では、ゲノムアセンブリは初めから染色体のコンセンサス配列として表現されている。大半のショートリードシーケンサーでは、ハプロタイプのバリエーションは捉えられない。
  5. セントロメアヘテロクロマチン領域の大半や、テロメア、低複雑度で長い繰り返し領域のシーケンスを決定することは、ヒトやマウスなど巨大プロジェクトでも困難である。
  6. ゲノムアセンブリにはエラーが存在し得る。例えばシーケンスエラーだったり、シーケンスブロックの順序の誤りなどが挙げられる。
  7. ゲノムアセンブリは一連のアセンブリヒューリスティックの結果として得られた結果にすぎないので、このデータをworking hypothesis(作業仮説 *1)にし、さらなる研究の基礎として活用すべきである。
  8. 1-7より、ゲノム配列が決定されるということは、当該の種のコンセンサスの配列が決定されることを意味する。

 

ゲノムシーケンシングとアセンブリの原則

  1. 現在、大半のゲノムプロジェクトはショットガンシーケンシングによって行われる。最初にゲノムDNAはランダムカットされ、所定のサイズ分布に断片化される。断片化したDNAからライブラリが作成され、シーケンスされる(論文図2)。
  2. 得られたシーケンシングデータは、適切な前処理を受け、洗練されたコンピュータアルゴリズムによってアセンブリされcontig配列になる。De novoアセンブリにはシーケンス断片間の十分なオーバーラップが必要であり、そのため、高いシーケンシングカバレッジ(またはリードデプスとも言う)でシーケンスが行われる。contig配列を繋げ、より長いscaffolfds配列が作成される。

 

ゲノムシーケンシング

  1. 最近(2014年)Finishした全ゲノムシークエンシングプロジェクトから判断すると、従来のサンガーシーケンシング(~1kb)とRocheの454シーケンシング(~800bp)からIllumina HiSeq(現在は通常150 bp)やSOLiD(通常50 bp)(SOLiDは既に販売停止)への高スループットシーケンシングへの移行が顕著である。
  2. Pacific Biosciences(最大5 kb)、IonTorrent(約500 bp)、Oxford Nanopore(論文では記載なし)など、より長いリードのシーケンシング技術が市場に参入している。
  3. 最近の研究では、異なるリード長のデータといくつかの異なるシーケンシングプラットフォームとの組み合わせによるアセンブリ方法が報告されている(Koren et al., 2012)。このようなハイブリッドアセンブリが常に単一のデータによるアプローチより優れているかどうかは依然として未確定だが、各方法の欠点が相殺されるため、この戦略は直感的により優れている。
  4. 従来のBACクローン単位でのサンガーシーケンシングでは、かなり限られたカバレッジでも大きな問題にはならなかった。しかし、ショートリード技術のみを使用する場合、高いリードカバレッジが必要になる。データが少なすぎるとアセンブリが断片化し、アノテーションやバリアントコールなどのダウンストリームアプリケーションで深刻な問題が発生する。
  5. アセンブリは、通常、大量のペアエンドのショートリードデータから始める(ロングリードだけのアセンブリ技術も存在する)。続いてコンティグをスキャホールドに融合させるためには、3〜40kbの長いインサートサイズを有する追加のライブラリーを生成する必要がある。どのような配列データが各ライブラリーのタイプおよびインサートサイズに必要であるかは、ゲノムのサイズおよびリピート含有量、ヘテロ接合性の程度およびアセンブリー品質(Sims et al、2014)を含む多くの要因にに依存するため、シーケンシングプロジェクトごとにパラメータを固有に決定する必要がある。
  6. 哺乳類ゲノムの大まかなガイドラインとして、45xカバレッジのペアエンドショートリード、45xカバレッジのミドルインサートサイズのライブラリ(3-10kb)、1-5×カバレッジのラージインサートサイズのライブラリ(10〜40kb)が提案されている(Nagarajan and Pop 2013 pubmed)。
  7. シーケンシングエラーの絶対数が増加するためカバレッジが高すぎることも問題になり得ることに注意する。著者ら自身の経験によれば、ショートリードのペアエンドライブラリの100xから50xへのカバレッジのダウンサンプリングは、アセンブリのいくつかのステップを大幅に改善しうる。
  8. 上記のパラメータを考慮するには、ゲノムサイズ、シーケンスエラー率、リピート含有量およびゲノム重複に関する基礎知識が必要になる。プロジェクトの開始時に目的の種の情報が得られない場合は、まずシングルエンドまたはショートインサートシークエンシングを使用して小規模なパイロットスタディを実施することを勧める。上記のパラメータは、k-merカウントアプローチを用いて近似することができる(Marçaisand Kingsford 2011; http://josephryan.github.com/estimate_genome_size.pl  Estimate genome size)。
  9. k-merカウントの実行および解釈方法に関する情報は、seqanswersなどのWebフォーラムで見つけることができる。一般に、ゲノムのリピート含有量が高い場合、正確なアセンブリのために、より長いインサートサイズのシーケンスデータが必要となる。ゲノムサイズの見積もりは、オンラインデータベースでも利用できる(box2 論文リンク

 

ウエットの作業

 ゲノムシーケンシングのウェットの作業はシークエンシングセンターに委託されることが多いので、この論文では、プロジェクトの計画段階で考慮すべき重要なライブラリー作成の基本的な手順と、下流の分析手順について、非常に簡単にだけ触れる。

  1. 個体のヘテロ接合性部位は、アセンブリに悪影響を及ぼす。高度に倍数性の種では、アセンブリは特に困難であり、特別に調整されたアセンブリパイプラインを必要とすることがある(Schatz et al、2012)。利用可能な場合は、近親交配の個体、単為生殖、または雌性発生後の子孫を使用することが一般的に推奨されている。個体の身元、年齢、性別、サンプリング時間および正確な場所など、将来の参照のために重要なメタデータでなければならない(Genome 10K Community of Scientists 2009 pubmed)。
  2. シーケンシングに使う組織: エネルギー的に活性な組織(筋肉など)は、高い割合のミトコンドリアDNA(mtDNA)が含まれ、アセンブリステップに問題を引き起こすリスクがあるため避けるべき。ターゲット生物以外の非標的DNAの度合いが高い腸や皮膚などの組織も避けることが推奨される。
  3. アセンブリ前にmtDNA配列リードを除去し、除去したリードはミトコンドリアゲノムのデータとして残しておく。それ自体が保存遺伝学にとって重要な情報を提供するかもしれない。ミトコンドリアゲノムを決めるには、除去したデータのうちほんの少しだけ使いアセンブリを行う。

 

DNAの品質

  1. 全長ゲノムシーケンシング、特にロングインサートサイズのライブラリーの場合、十分な量の高品質で無傷の非分解DNAが必要になる(Wong et al、2012)。複数ライブラリーを使用する全ゲノム解析では、starting materialとして約1mgのDNAが必要になる(ショートインサートのペアエンドライブラリの場合約6μg、2-10kbライブラリーの場合約40μg、> 20kbライブラリーの場合約60μg: 論文執筆時点の話) 。
  2. 高品質のDNAを大量に取得することが不可欠であるが、これは保全上の懸念がある多くの種にとって大きな障害となりうる。捕獲可能な動物が利用可能である場合、しばしば高品質DNAの供給源として利用できるが、そのような供給源から同定されたゲノム変異は野生集団を代表するものではないことに注意する。 DNAサンプルを提出する前に、高解像度ゲル(例えばパルスフィールド電気泳動;試料は典型的に> 100kbの断片を示すはずである)上でその完全性を調べるべきである。

 

ライブラリの準備

  1. 現在のライブラリ調整のほとんどの技術は無視できない回数のPCRに依存した方法のため、無視できない数のduplicartion readsが発生することに注意する。ただし、カバレッジが非常に高い場合、偶然全く同じ部位を読んだリードが二重に発生する可能性があり、これはPCR duplicartionと区別できない(*2)。
  2. duplicartion readsには付加価値がなく、カバレッジベースのクオリティチェックを損なう可能性があるため、アセンブリ前に削除する必要がある。重複は一般的に、ショートインサートサイズのライブラリー(<500bp)の数パーセントを構成する。ロングインサートサイズのライブラリー(> 10kb)では95%を超えることがある。
  3. 別の課題は、ペアエンドライブラリのインサートサイズの決定である。一般的に0.2〜40kbの範囲の大きさで複数使い、短いインサートサイズのライブラリーのカバレッジが多いことが望ましい(Gnerre et al、2011)。 20kbを超えるインサートサイズは、アセンブリの最終的な連続性に大きな差異をもたらすが、高品質で生産するのは容易ではなく、現在多くのシーケンシングセンターで制限因子となっている。
  4. もう一つの重要な問題はペアエンドシーケンスのリードの向きである。使用される技術と元のDNA断片に関連して、リードは内向き(→ ←;例えばイルミナペアエンドシーケンシング)または外向き(← →;例えばイルミナメメアイトペアシーケンシング)がある。向きが異常なペリードは、予想外に短いインサートサイズ由来する可能性がある。また、異常な向きやサイズのメイトペアは、隣接していないゲノム領域のキメラであることが多い。
  5. このようなアーチファクトデータは、アセンブリ前にフィリタリングで除く必要がある、しばしば、トリミング後に組み立てるために使用可能な独特の読み取り対のほんの一部しか残さない。データを正しく処理するためには、データを扱う際に常にライブラリのことを分かった上で('library aware')作業する必要がある。

 

データのマネージメント(重要 *3)

  1. 通常のゲノムシーケンシングプロジェクトで生成されるデータの量は驚異的である。 100×カバレッジを有する脊椎動物のゲノムでは、数百ギガバイトのオーダーのデータファイルになる。アセンブリ工程中、テンポラリファイルはテラバイト境界を容易に越える。よってプロジェクトの開始時には、すでに適切なデータ管理とバックアップ戦略が取られていなければならない。多くの大学は、データ保管施設を含む地方または国のコンピューティンググリッドに接続されており、可能であればこれらを活用することが強く推奨される。コンピューティングインフラストラクチャで働くバイオインフォマティクスの専門家は、生物学の研究者とコンピューティンググリッドシステムの専門家との間の重要なリンクを提供する(Lampa et al。2013 pubmed)。このようなコラボレーションは、プロジェクトの計画段階で既に確立されているはずである。
  2. シーケンシングセンターは、多くの場合、データ分析とアセンブリの支援も提供する。しかし、それらの自動化パイプラインは、非モデル生物のデータには最適化されず、保存生物学の観点からは有用ではない可能性がある。したがって、プロジェクトの開始前にどのような支援が施設によって提供されるかを明示的に議論することは極めて重要である。より一般的には、アセンブリの計算ステップを実行するのに十分な専門知識がコア研究グループに存在するかどうかを検討する必要がある。大規模シーケンシングデータの大部分のデータ処理およびゲノム解析は、UNIXベースのオペレーティングシステムを実行する高性能コンピューティングクラスタ上で実行される。全ゲノムシーケンシングデータを処理するには、バイオインフォマティクスの専門家である必要はないが、UNIX環境とコマンドラインソフトウェアに関する十分な基本知識、シェルスクリプトの作成、生物学的データ分析のための一般的なスクリプト言語PerlPythonなど)を適用できる必要がある。

 

アセンブリ前の作業

  1. アセンブリに先立って、シーケンシングデータの品質、GC含量、リピート存在量または重複した読み取りの割合を評価する必要がある。要約統計を提供するFastQC(http://www.bioinformatics.babraham.ac.uk/projects/fastqc)のようなツールは、有用な出発点である。
  2. PCR dplicationによる低品質データのトリミングは、様々な異なるソフトウェアおよびスクリプト(例えば、ConDeTri; Smeds andKünstner 2011)を用いて行うことができる。 k-merカウントアプローチ(例えばSOAPdenovoパイプラインで適用される)を使用したスタンドアロンのエラー修正は、多くのデータセットの代替として有用である。ただし、品質フィルタリングの最適な厳密性は、個々のプロジェクトと対象となるアセンブリパイプラインに固有のものであることに注意する。
  3. アセンブリパイプライン内でトリミングとエラー訂正が実行されるALLPATHS-LG(Gnerre et al、2011)のような一部のアセンブラでは、トリミングしたデータでなくオリジナルのシーケンシングデータを使う必要がある。
  4. ライブラリーのプライマーおよびベクター配列は(シーケンシングマシンがそれらを除去したと主張したとしても)データ中に残存する可能性が高い。単純なスクリプト(cutadapt; Martin 2011)で除去することができる。また、Illuminaシークエンシングでは、シーケンスクオリティスコア較正のために、PhiXファージ由来のDNAをシークエンシングに加えることが多い。そのような豊富な汚染配列を除去しなければ、(核ゲノムと比較して高いカバレッジのために)アセンブリプロセスが中断され、キメラおよびコンタミネーションのコンティグが生じる可能性がある。
  5. raw データから既知のベクター汚染を除去する最も簡単な方法は、ショートリードアライナー(BWA; Li and Durbin 2009)を使用し、汚染配列にマッピングされる全てのリードを捨てることである。

 

 De novo assembly

  1. ゲノムアセンブリのツールは様々存在するが、速度、スケーラビリティ、最終配列の品質(Miller et al、2010; Earl et al、2011; Narzisi and Mishra 2011; Bradnam et al、2013)の性能において大きく異なっている。アセンブリ方法の中には他の方法よりも明らかに優れているものもあるが、特定の状況でどのツールが最も適切かを予測することは現在困難である。
  2. すべてのアセンブリプロジェクトは、生成されたデータ構造、サイズ、塩基組成、リピート含有量、多型のレベルなどの点で固有である。
  3. ショットガンシーケンシングデータのデノボアセンブリに利用可能なソフトウェアが多数あり、常に新しいプログラムがリストに追加されている。いくつかのアルゴリズムはミス・アセンブリを最小限に抑えることに焦点を当てているが、他のアルゴリズムは主に連続性(contiguityのこと)を改善することを目的としている(時には精度を犠牲にして)。
  4. ほとんどのアセンブリアルゴリズムは、ライブラリサイズの分布に従って最適に実行されるため、プロジェクト計画やシーケンシングの段階で既にアセンブリ戦略の選択を検討することが重要になる。主な文献やアセンブリソフトウェアのウェブサイトの情報、さまざまなウェブフォーラムが、最新のディスカッションや他の研究者の経験を共有するための入り口として役立つ(box2 論文リンク
  5. 従来のSangerシークエンシング(例えばCeleraアセンブラ、ArachneおよびPCAP; Batzoglouら2002; Huangら2003; Denisov et al、2008)またはRoche 454シークエンシング(例えば、 Newbler)は、オーバーラップレイアウトコンセンサス(OLC)と呼ばれるアセンブリアプローチを使用している。OLCアルゴリズムは、IlluminaまたはSOLiDデータでは、一般的に計算量が多すぎるとみなされる(too computationally intensive)(時間的に不利)。それでも、Edena(Hernandez et al、2008)、SGA(Simpson and Durbin、2012)、FERMI(Heng Li、2012)などのいくつかのアセンブラはOLC戦略を追求している。
  6. ショートリードのデノボアセンブリのほとんどの戦略は、extension-basedの方法とDe Bruijn(またはオイラーの)グラフアルゴリズムの2つのクラスに大別することができる(Nagarajan and Pop、2013)。extension-based methodsのアセンブラの SSAKE(Warren et al、2007)やJR-Assembler(Chu et al、2013)はメモリと計算時間両方で計算上非常に効率的だが、シーケンシングエラー、反復領域および高レベルのヌクレオチド多型(Chu et al、2013)に非常に敏感に影響される。
  7. 現在、最も一般的に使用されるアプローチはDe Bruijnグラフに基づいている。リードは長さkのシーケンリースリードの部分文字列k-mersに変換され、 k-1 merを共有するときにリンクされグラフ(ネットワーク)のノードを形成する。 SOAPdenovo(Luo et al、2012)、ALLPATHS-LG(Gnerre et al、2011)、ABySS(Simpson et al、2009)、Velvet(Zerbino and Birney 2008)などの高度に使用されたアセンブリソフトウェアは、すべてDe Bruijnグラフアルゴリズムを採用している。
  8. 異なるアルゴリズムの特徴を組み合わせ、複数のシーケンシング技術のデータを利用する ‘hybrid’なアルゴリズムも存在する: Atlas(Havlak et al、2004)、Ray(Boisvert et al、2010)、MaSuRCA(Zimin et al、2013) 。
  9. 一般的には、いくつかの異なるアセンブリ方法をテストし、現在の特定のデータに最も適しているものを評価することが推奨される。
  10. ドラフトゲノム構築の過程では、数回のアセンブリ、評価、パラメータ調整が繰り返されるiterativeなプロセスとして扱われなければならない。さまざまなアセンブリアルゴリズムとソフトウェアのより包括的なレビューについては、例えば、Miller et al、2010 pubmed; Nagarajan and Pop 2013 pubmed を参照。
  11. 最初のコンティグ構築後、コンティグをScaffoldsに結合するために、ロングインサート(mate-pair、fosmid-end、jumping)ライブラリー(Zhang et al、2012)からのペアリード情報を用いるのが一般的である。
  12. 例えば短い低複雑度の領域(low-complexity regions)をブリッジするために、追加のショートインサートペアエンドライブラリはしばしば有用である。
  13. コンティグ間のギャップの長さは予想されるインサートサイズから推定され、通常Nsのストレッチで満たされる。Scaffoldingステップは、一般的に使用される多くのアセンブリプログラムにすでに含まれているが、SSPACE(Boetzer et al、2011)やBESST(Nystedt et al、2013)などの独立したアプリケーションもある。
  14. Scaffoldingで生じたNsは、元のペアエンドシーケンシング情報を使用して、GapCloser(Li et al、2010)、GapFiller(Boetzer and Pirovano、2012)、iMAGE(Tsai et al、2010)などのソフトウェアで除去できる。最近では、例えばPac​​Bioなどのロングリードシーケンシングデータを使う方法も出現している(English et al、2012)。
  15. アセンブリソフトウェアを選択する際には、シーケンシングデータの量と使用可能な計算リソースの両方を考慮することが重要になる(Schatz et al、2010a)。 SOAPdenovoやALLPATHS-LGなどのDe Bruijnグラフメソッドは、一般に大量の計算用メモリ(RAM)を必要とする。哺乳類サイズのゲノム(〜3Gb)のアセンブリでは、シーケンシングデータの量に応じてテラバイトの内部メモリが必要となる(Lampa et al、2013)。大規模なコンピュータクラスタがローカルで利用できない場合は、共同機器の購入、バイオインフォマティックスグループとの共同プロジェクト、または市販のコンピューティングクラウドの利用を検討する必要がある(Box 2; Schatz et al。2010b pubmed)。
  16. もう一つ考慮すべき点は、自由に利用できるプログラム(このカテゴリーに含まれるほとんどのプログラム)を使用するのか、商用ソフトウェア(CLCワークベンチやDNASTARのLasergeneなど)に投資するのかである。市販のソフトウェアは、通常、自由に利用できるプログラムよりもユーザーフレンドリーであるため、限られたバイオインフォマティクススキルを持つ研究者も容易に使用できる。商用ソフトウェアの欠点は、購入やライセンス供与に伴う(しばしば実質的な)コストを除けば、アルゴリズムの詳細を調べたり変更したりすることは不可能に近い「ブラックボックス」ソリューションに似ていることである。一般に使用されるソフトウェアアプリケーションの中には、シーケンシング機器と一緒に配布されるものもある。
  17. Highly closedな生物分類群のゲノム情報が利用できる可能性が存在する(Kohn et al。2006)。最も一般的な方法は、最初にコンティグを新たに作製し、次いでそれらをreferenceのゲノムにアライメントさせscaffolding を助けることである。広範なシンテニーおよび遺伝子の順序の保存を仮定すると、そのようなアプローチは、低いカバレッジデータ(Kimら、2013)や非常に短いショートリードシーケンシング(Wang et al。2014)を用いても大きなscaffoldsを構築することを可能にする。別のアプローチは、いわゆるAlign-Layout-Consensusアルゴリズムである。この方法では、デノボアセンブリのオーバーラップ工程は、リードのhighly closedなリファレンスゲノムへのアラインメントに置き換えられる。次いで、互いにオーバーラップするリードの情報を使用して新たにコンティグおよびscaffoldsが構築される(Schneeberger et al、2011)。

 

クオリティ評価

  1. アセンブリが正常に終わった後も、ユーザーは品質を評価したり、さまざまな方法で複数のアセンブリを比較したりすることが必要になる。しかし、上で論じたように、すべてのドラフトゲノムアセンブリは真のゲノム配列の仮説を構成するだけであり、真実を知ることができなければ、その質の評価方法には課題が残る。
  2. アセンブリのさまざまな側面を反映するさまざまな指標が利用可能である(Bradnam et al、2013 pubmed)。それらは外部データからの追加情報を必要とするアプローチと、アセンブリ自体の情報に基づくものだけに大別される。
  3. 1つの基本的なメトリックは、アセンブリに含まれるゲノムの割合である。予想されるゲノムサイズは、C-valuewiki、他参考HP)、またはk-merベースのアプローチから推測できる。アセンブリの連続性を評価するためのもう1つの標準的な測定基準は、N50 statisticである:定義上、組み立てられたヌクレオチドの50%は少なくともこの長さのコンティグ(コンティグN50)またはscaffolds(scaffolds N50)に見出される。したがって、N50 statisticはアセンブリされた配列の一種の中央値を記述する。
  4. N50のバリエーションとして、最近では期待されるゲノムサイズを組み込んだNG50やNG Graph(Bradnam et al、2013)が導入され、異なるアセンブリ間のcontiguity を視覚化して評価する効率的な手段が提供されている。
  5. しかしN50 statisticとその関連statisticは単に連続性を示し、アセンブリの精度に関する情報は含まれないため慎重に解釈する必要がある。アセンブリ内のエラーを検出するために、ペアエンドまたはペアペアデータを、再マッピングしてその情報を使用することができる(例えば、ソフトウェアREAPR; Hunt et al、2013)(紹介)。カバレッジが低い領域または誤った向きでペアエンドリードがマッピングされればミスアセンブリが示唆され、異常なインサートサイズは小さな挿入または欠失がアセンブリに起きていることを示す。
  6. その領域だけ非常にカバレッジが高く、SNPが非常に多かったり、また大部分のリードで同じ塩基だが他の複数リードで別の塩基を表すような領域は、ほぼ同一の(だが崩壊した)リピートの存在を示している。
  7. これらを検討するソフトウェアアプリケーションは数多くあり、その例は現在の文献(Earl et al、2011 pubmed; Bradnam et al、2013 pubmed)に記載されている。The amosvalidate pipeline(Phillippy et al、2008; Schatz et al、2013)は、1つのパイプラインでいくつかのゲノムアセンブリ診断を包括しているがsmall~middleサイズのゲノムに最適化されている。
  8. 該当生物に関して独立した実験データセットがあれば、おそらく、外部情報の最良の供給源となる。例えば、optical maps(wikiショートリードの憂鬱)からのデータは、scaffoldsの精度を検証し、scaffoldsをさらに拡大して染色体レベルに近づけることを可能にする。同様に、BACまたはfosmidライブラリーを個々にアセンブルされた配列は、配列の正確性およびリピート含量を評価するのに役立ち得る。しかし、どちらのアプローチも、正確なアセンブリ自体に依存しており、現在小規模な研究所では容易に入手できない。
  9. RNA抽出のための新鮮な組織にアクセスすることは深刻な制限かもしれないが、ショットガントランスクリプトームシーケンシングデータ(RNA-seq)からの独立したデノボアセンブリはより容易に生成でき、ESTはすでに実施されているかもしれない。したがって、transcriptome データの内容およびexon構造は、遺伝子がscaffoldsをまたいでいる場合、scaffoldsの修正のために重要な追加リソースとなり得る。
  10. 比較ゲノムのアプローチは、追加のデータの生成を必要としない別の手段を提供する。例えば、オーソロガスなcore eukaryotic protein sequencesの存在および完全性の(Parra et al。2007)評価は、アセンブリが包括的であるかに関する最初の直感的情報を提供してくれる。
  11. Sister taxaの高品質のリファレンスゲノムが存在するなら、ゲノム比較することは、ブロードなレベルでのシンテニーおよび遺伝子順序の保存を仮定して、ミスアセンブリおよびキメラコンティグを検出するガイダンスになり得る。しかし、小規模なrearrangementsは現実的であり、深い調査を必要とするかもしれない。
  12. 他の生物由来のDNAが、様々な段階で(サンプリングおよび実験手順の両方の間に)ゲノム試料を汚染している可能性がある。サンプリング段階での汚染は、研究種に関連する寄生虫または他の微生物に関する情報を運ぶ可能性もあり、実際には保存の観点から興味深いかもしれない。
  13. 外部ゲノムリソースがあればそのようなコンタミネーションを見つけるのを助けてくれる。このようなコンタミネーションの痕跡を見つけるために、BLAST検索または類似のローカルアライメントがしばしば用いられるが、結果は注意深く解釈される必要がある。正確にアセンブリされた配列であっても、特にターゲット分類群内のシーケンシングが不十分である場合、よくアノテーションづけ付られたゲノムを有する遠い種からベストヒットが出る可能性がある。同様に、scaffoldsの大部分がターゲットクレードに顕著なヒットをもたらしていると、scaffolds内の小さなストレッチでのコンタミネーションが完全に見逃される可能性がある。
  14. ヒトゲノムが他の哺乳動物ゲノム配列に見つかる場合、サンプルの取り扱いに原因があると解釈するのが一般的であり、特に問題となる。しかし、哺乳類ゲノムの一部については、該当領域の配列が決定されていないため、シーケンシング精度がもっとも高いヒトまたはマウス配列に対してBLASTのベストヒットが出ているだけかもしれない。

 

ゲノムアノテーション

  1. ゲノム配列のポテンシャルを完全に引き出すには、GOterm(遺伝子オントロジー)(Gene Ontology Consortium 2004; Primmer et al、2013)、'Kyoto encyclopedia of genes and genomes' (KEGG) (KanehisaおよびGoto、2000)、microRNA and epigenetic modifications (ENCODEプロジェクトコンソーシアム2012)のような遺伝子モデルおよび機能情報からなる生物学的に関連する情報をつけてアノテーションする必要がある。
  2. 非モデル生物では、アノテーションはしばしばタンパク質コード配列(CDS)またはtranscriptsに限定されることがより一般的である。既存の遺伝子モデルがほとんど欠如しているため、新たにシーケンスされた種に遺伝子にアノテーションを付けるにはかなり難題だが、自動アノテーションは個々の研究グループにとって可能になった(Yandell and Ence 2012)。それでも、完全なゲノム注釈はかなりの労力を要し、バイオインフォマティクスの熟練を必要とする。一般的なワークフローについてのみ説明する。包括的なレビューはYandellとEnce(2012)を参照。
  3. アノテーションが成功するには、ゲノムアセンブリの品質に強く依存する。小さなギャップによってのみ切断された、ほぼ完全なゲノム(~90%)のみ満足のいく結果が得られる。経験則として、大きなゲノムはより長い遺伝子を有しており、したがって、アノテーションを成功させるためには、より連続したアセンブリが必要である(YandellおよびEnce 2012の論文の図1参照)(pubmed)。
  4. アノテーションプロセスは概念的に2つのフェーズに分けることができる。「計算フェーズ」は、他のゲノムからの複数行の証拠または種特異的なトランスクリプトームデータを並行して使い、初期の遺伝子とtranscriptsを予測する。第2の「アノテーションフェーズ」では、アノテーションパイプラインによって決定されたルールのセットに従って、すべての(時には矛盾する)情報が遺伝子アノテーションに合成される。
  5. 遺伝子予測の前に、複雑性の低い領域および転移因子を含むリピート配列をマスクすることが極めて重要である。リピートは種間で保存されていないことが多いため、RepeatModelerやRepeatExplorerなどのツールを使用して種特異的リピートライブラリを作成することが推奨される(Nováket al、2013)。
  6. リピートがマスクされると(例えば、RepeatMasker; http://www.repeatmasker.org)、関連種の遺伝子モデルを使いトレーニングしたab initioなアルゴリズムをコード配列(CDS)のベースライン予測に使用することができる(例えばAUGUSTUS; Stanke et al、 2006)。
  7. タンパク質アラインメント(例えば、tblastxを用いる)および種々の他の種からの合成タンパク質リフトオーバーは、遺伝子モデルの予測を補う貴重なリソースを提供する。
  8. 最良の証拠はESTとはRNA-seqデータから得られる。これらから、CDSに加えて、スプライスサイト、転写開始部位および非翻訳領域(UTR)に関する情報を遺伝子モデルに提供される。可能であれば、mRNAはストランド特異的にシーケンシングされる必要がある。
  9. 第二段階では、最初の予測およびタンパク質、ESTまたはRNAアライメントからのすべての証拠を、遺伝子アノテーションの最終セットに合成する。証拠はほとんど不完全であり、また時には矛盾するので、これは手作業のキュレーションなしには難しい作業になる。それでも、MAKER(Cantarel et al。2008)やPASA(Haas et al。2003)のような自動注釈ツールがいくつか存在し、そこからの証拠をウエイトをつけて取り入れ評価している。
  10. 9のようなツールは一般に良好な結果をもたらすが、定性的検証は重要である(例えば、オープンリーディングフレームの長さを評価することにより)。アノテーションの目視による検査は、イントロン漏れ(intron leakage)(イントロンがpre-mRNAの存在のためエキソンとして注釈される)または遺伝子融合などの系統的な問題を検出するための別の必須のコンポーネントである。
  11. GMODプロジェクトのWebApollo(Lee et al、2013)(ユーザーガイド)のようなツールは、ユーザーがビジュアルインターフェイスを介して直接アノテーションを編集できるようにするため特に便利である。

 

ゲノムの公開

  • ドラフトゲノム配列は現在、ますます増加している。伝統的なEMBL(European Molecular Biology Labs)のENSEMBLやWellcome Trust Sanger Institute、National Center for Biotechnology Information (NCBI) のデータベースではゲノムやメタ情報へのアクセスを提供するが、既に全ての配列のアノテーションやキュレーションは出来なくなっている。
  •  したがって、NCBIはすでに、ユーザが生成したゲノム配列とアノテーションのドラフトをアップロードを受け付けている。 他のユーザーがアセンブリとそのアノテーションを改善できるようにするには、すべてのrawデータをNCBIのBioProjectなどにアップロードする必要がある。

 

Perspectives

 ゲノミクス分野に進出した保全または生態学研究グループの実践的指針を提供する目的で、全ゲノムシーケンシング、アセンブリアノテーションの現在の方法に関する情報を要約した。 焦点は、保全の観点から、非モデル生物の大きくて複雑なゲノムに焦点が当てられている。 はじめに、一般的なゲノムリソースと、完全なゲノム配列が保存生物学の設定に適用できるいくつかの異なる方法を概説した(論文 図1も参照)。 保全ゲノム学は若い分野であり、適用された保全状況においてゲノム資源が試験に供されている例は依然として限られているが、そのようなケースのいくつかを述べる。

  1. イルミナのシーケンシングテクノロジーを用いられた最初の非モデル生物ゲノムの1つはジャイアントパンダ(Li et al。2010)であった。パンダゲノム論文の焦点は保存に関する問題ではないが、フォローアップ研究では、ゲノムドラフトを利用して、人口構造、適応的な遺伝的変異および人口統計を推測した(Wei et al、2012)。
  2. 同様に、アイマイでは、マダガスカルの異なる地域の12人からのリシークエンシングデータを利用して詳細な遺伝子集団構造を推定し、 landscape genetic analysesを実施した。この結果は、マダガスカル北部の大規模かつ連続的な生息地を維持するための保全資源の配分に関するガイダンスを提供するために用いられた(Perry et al、2013)。ゲノム資源は、伝染性の顔面のガンのため、野生で絶滅の危機に瀕しているタスマニアデビルの育種プログラムに利用されている。リファレンスゲノム配列とゲノムワイドリシークエンシングデータの組み合わせにより、腫瘍形成に関与する候補遺伝子の同定を含む、この疾患の多くの詳細を調査することが可能になった(Murchison et al、2012)。
  3. カリフォルニア州のコンドルの繁殖プログラムにおける突然変異を引き起こす発生性疾患の広がりを制限するために利用された(Romanov et al、2009)。最終的に、ゲノムワイドSNPスクリーニングは水産資源のモニタリングおよび管理に関するいくつかの研究において有効であった(Primmer 2009; Nielsen et al、2012a)。

 

今後の方向性

 ナノテクノロジーのシーケンシングの急速な進歩と計算方法のさらなる発展により、ワークフローのすべてのステップが引き続き改善されると期待される。新しいライブラリー調製プロトコルにより、より少ないstarting materialからのシーケンシングが可能になり、より正確かつより正確に推定されたインサートサイズを有するライブラリーが生成され、エラー率が低下してより長いリードが生成される。より効率的なアセンブリアルゴリズムの開発と計算能の向上は、バイオインフォマティクスのデータ処理をより広い範囲のユーザーに受け入れやすくするだろう。ゲノムシーケンシングとアセンブリに伴うコストが低下し続ける中、ドラフトゲノム配列の生成は、ゲノムの大きな種についてもまもなく日常的になるだろう。この発展により、限られた資金しか持たない小規模な研究グループでさえも、保全生物学および関連する分野におけるゲノムアプローチの使用がエンハンスされ、選択した種のゲノム資源の発展が可能になる。

 進歩のもう1つの重要な領域は、低品質な博物館資料から得られたサンプルの非侵襲的サンプリングによる、時間を超えたゲノム解析の発展である。ゲノムデータを保存し共有する方法を開発することも、これらの資源を保存のために重要になる。これらの有望な発展にもかかわらず、科学だけでは将来の保全の課題を満たすには不十分であることを認識する必要がある。したがって、保存遺伝学からゲノムスケールデータへの技術的移行は、適用される保存生物学がどのようにゲノムデータから最も利益を得ることができるかについての議論を緊密に伴う必要がある(McMahon et al、2014 pubmed)。この議論は一般的なケースバイケースで行われる必要があり、科学者や政治意思決定者も参加する必要がある。

 

引用

A field guide to whole‐genome sequencing, assembly and annotation

Ekblom R, Wolf JB.

Evol Appl. 2014 Nov;7(9):1026-42

 

 参考

*1

サイエンスの大前提は反証可能性があること。

 

*2

Amplicon seqやRNA seqでは偶然の完全一致が無視できない数で起こりえるため、PCR duplicationを除くのは難しい。

 

*3

このブログでは、この段落の重要性を特に強調したい。

高速な端末エミュレータ Alacritty

 

  AlacrittyはRustで書かれたGPUレンダリングに使う高速な端末エミュレータOpenGLwiki)を使ってレンダリングを行う。開発はまだアルファ段階らしいが、すでに色々なプラットフォームに対応している(windowsはこれかららしい)。開発の大きな動機は、WUXGAから4k、5kのような高解像度ディスプレイでtmuxを起動し、vimを開いてコーディングしている時の画面描画の遅さだったと述べている。

 

Announcing Alacritty, a GPU-accelerated terminal emulator

https://jwilm.io/blog/announcing-alacritty/

 

Alacrittyに関するツイート。

twitter.com

 

インストール 

mac os10.13に導入した。

本体 Github

mac osへのインストール

#Rust compiler導入(画面の支持に従い導入する。既に入れているなら以下の三行不要。)
curl https://sh.rustup.rs -sSf | sh
export PATH="$HOME/.cargo/bin:$PATH" #<=をbash_profilenに追記
source ~/.bash_profile #設定を反映

#本体導入
git clone https://github.com/jwilm/alacritty.git
cd alacritty
#ビルド
make app
cp -r target/release/osx/Alacritty.app /Applications/

または参考ページ(*1)のようにbrewを使って導入する。/Applications/Alacritty.appをダブルクリックして起動する。または他の端末エミュレータ環境から呼び出す。

 

使用例

1、準備

mac pro2012モデル(Nvidia GTX680 搭載)でテストした。画面のレンダリング負荷が高くなる状況を作り出し、Alacrittyとターミナルで作業に影響するほど差が出るのが調べてみる。ディスプレイは4k一枚だけとする。

Nvidiaは昔からopen GLに弱いと言われている。はじめにCinebenchレンダリングによるCPUとGPUベンチマークツール)でopen GLスコアをチェックしておく(*12)。OpenGLスコアが低すぎれば、そもそもAlacrittyの高速レンダリングは期待できないことになる。

f:id:kazumaxneo:20180815103223j:plain

Open GLのスコアは50.4だった。2018年現在では高いスコアでは全然ないが、一昔前のオンボードグラフィックほど極端に悪いわけでもない。このグラフィックカードでテストしてみる。 

 

-- 補足 --

ターミナルでコマンド入力時、カーソルの移動速度が遅いなら以下を最大にしてください。

f:id:kazumaxneo:20180815112046j:plain

文字の誤入力が増えたら微調整してください。その他、Ctrl+e(行末に移動)、Ctrl+a(行頭に移動)等 (*2) のキーバインドを使う癖も忘れずに。

 

2、Alacritty起動

/Applications/Alacritty.appを起動する。

f:id:kazumaxneo:20180814233144p:plain"bash --login"が自動実行されなかったので、tmuxなどにパスが通っていない。Alacritty初回起動時に自動作成される$HOME/.config/alacritty/alacritty.yml に以下を追記する(*3)(カラースキームの変更等もalacritty.ymlをいじる)。

bashの場合

shell:
program: /bin/bash
args:
  - --login

zshの場合

shell:
program: /usr/local/bin/zsh
args:
  - --login

 

3、tmuxと監視ツールの起動

Alacrittyを再起動し、tmux等のパスが通っているか確認する(tmuxインストールは*4参照)。O.Kならtmuxを起動する。

> tmux

f:id:kazumaxneo:20180815114903j:plain

tmux起動中は下に帯が付く。

 

横2つに画面分割する(ペイン分割と呼ぶ)。 Ctrl b、続いてshift %と打つ。

f:id:kazumaxneo:20180815115029j:plain

分割された。両方独立して動作する。左ペインがメインのコマンド実行パネルで、右ペインでファイル名を確認しながら進めることもできる。

 

ペインを指定して切り替えてみる。Ctrl b、続いてqと打つ(*5参照)。画面に番号が表示されるので、番号が表示されている間に右側ペインの番号を打つと、1が割り当てられている右側ペインに移る(左は0)。

f:id:kazumaxneo:20180815121144j:plain

カーソルが右側ペインに移っている。

 

ペイン1はhtopの監視画面とする(インストールは*4参照)。htopを起動する。

> htop 

f:id:kazumaxneo:20180815120223j:plain

htopはプロセスがモニタリングできるtopの高機能版(*7)。rootなら指定プロセスを選んですぐにkillできる(*8)。マウス操作もある程度受け付ける。

 

次はペイン1を縦2つに分割する。 Ctrl b、続いてshift "と打つ。

f:id:kazumaxneo:20180815121230j:plain

ペイン1が縦2つに分割された。ここでは右上がそのままペイン1、右下がペイン2になる。

 

Ctrl b、続いてq、続いてを押し、ペイン2に切り替える。ペイン2は vtopを起動してCPUとメモリのの監視画面とする(インストールは*4参照。解説は*9-10参照)。

> vtop

f:id:kazumaxneo:20180815121312j:plain

htopと似ているが、トータルのメモリ使用率とCPU使用率はこちらの方がわかりやすい。カラースキームはdefaultのままにしている。

 

キーボードのqを押すと終了する。更新頻度をかなり早めにしてみる。テーマも変える。例えば20ms/secに更新頻度を上げるなら

> vtop --theme wizard --update-interval 20s

f:id:kazumaxneo:20180815131427j:plain

CPUのdotの波の流れが速くなった。GPUを使って画面レンダリングしているので、動きが非常に滑らかになっている。

 

3.5、おまけ。デタッチとアタッチ

tmuxは本来、複数の計算機サーバを並行して使うときなどに活躍する(*6参照)。一度デタッチしてみる。Ctrl-b、続いて dを押す。

f:id:kazumaxneo:20180815130933j:plain

tmux起動前に戻っている。

 

tmux aと打てば直前にデタッチしたセッションに再度アタッチできる。

> tmux a

f:id:kazumaxneo:20180815131427j:plain

デタッチ前の状態に復帰している。複数セッションを使っている時はtmux list-sessions を押してセッション番号を確認後、" tmux a -t #番号"で指定セッションに復帰する。tmux起動中にCtrl-b、続いて wを押せば、プレビューを見て画面を切り替えることもできる。デタッチしてもセッションはバックグラウンドで動いており、例えばマッピングなどを実行中にデタッチしても、マッピングのジョブは動き続けている。

 

 

4、mac純正ターミナルとAlacrittyの比較

Ctrl b、続いてq、続いて0を押し、左のペイン0に切り替える。nvimでpythonのコードを表示してみる(インストールは*4参照)。

> nvim test.py

f:id:kazumaxneo:20180815134218j:plain

これで1ウィンドウに3つのペインが表示されたことになる。

 

mac純正ターミナルでも同じことを行い、Alacrittyと比較してみる。

f:id:kazumaxneo:20180815135232j:plain

右側が純正ターミナル、左側がAlacritty。

 

nvimで開いたコードのスクロール速度。右側のターミナルはテーマ "ICEBERG"(リンク)を入れてある(*11)。

キャプチャだとフレーム不足で分かりにくいが、右側の純正ターミナルではスクロールがカクカクしていて、レンダリングがキーの打ち込みに全然追いついていない。対して左側のAlacrittyでは、キーの打ち込みスピードにダイレクトに応答する感じで非常にスムーズにスクロールされる。主観では3倍くらいはスムーズに動いているように感じる。

 シロイヌナズナゲノムを表示する。

やはりキャプチャだと分かりにくいが、Alacrittyの方がずっとスムーズに動いている。100 fps以上に対応したディスプレイならもっとスムーズに動くかもしれない。

 

上に下に行ったり来たりしていて、スクロールルカクカクになるとストレスが溜まります。高解像度ディスプレイとvimなどの組み合わせて開発されている方は恩恵が大きいのではないでしょうか?どなたか忘れてしまいましたが、ツイッターで教えてもらいました(ホタペンさんでした?)。ありがとうございました。

引用

https://github.com/jwilm/alacritty

 参考ページ

*1  alacritty参考ページ

https://qiita.com/zebult/items/0047b72916383a5c0acf

*2

ターミナルでのショートカットキー - 総合的な学習のお時間

*3

ログインシェルとインタラクティブシェルと~/.bashrc達の関係

*4 今回使用したソフトのワンライナーインストール

#brewリンク
/usr/bin/ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"

#tmux
brew install tmux

#vtop
brew install npm #Node.jsのパッケージマネージャ。持ってない人は先に導入。
npm install -g vtop

#htop
brew install htop

#neovim
brew install neovim

*5

tmux起動時に画面分割してサイズ調整させる

*6 tmuxの本来の使い方

インフラエンジニアならtmuxを使いこなしているか!? - Goalist Developers Blog

*7 htop解説

https://qiita.com/khotta/items/a998a52bc56c3bfe5f91

*8 htop操作 

https://www.submit.ne.jp/379

*9 vtop公式

https://parall.ax/blog/view/3071/introducing-vtop-a-terminal-activity-monitor-in-node-js

*10 vtop解説

https://orebibou.com/2017/05/ターミナルでグラフィカルに監視する『vtop』/

*11 iceberg

気分転換にTerminal.appの配色を変えてみる?暗青系のテーマ「Iceberg」を移植しました - ここぽんのーと

*12

Cinebenchのテストは、他の環境と比較できるよう固定解像度(=固定ピクセル数)で実行される。ディスプレイの見かけの画面解像度を変えてもベンチマークスコアは変化しない。

 

そのほか

Macで使えるターミナルエミュレーターアプリ達 

https://rcmdnk.com/blog/2017/02/14/computer-mac/

 

アライメントビューア Alan

 

Alanはターミナル(端末エミュレータ)で動くアライメントのビューア。GUIに頼らず端末内で確認作業を完結させることができる。

 

インストール

ubuntu18.04でテストした。

--対応フォーマット--

FASTA、Clustal format alignments

本体 GIthub

alan-2.1.1/alangを叩いて実行する。bashスクリプトになっている。

 ラン

マルチFASTAファイルを開く。入力フォーマットは自動認識する。明示する時は"-n"(DNA)か"-p"(プロテイン)をつける。

alan input.fa

配列が端末エミュレータ(ターミナル)内に表示される。画面からはみ出している場合、十字キーで上下左右に移動できる。

f:id:kazumaxneo:20180814120119j:plain

 

Mafft紹介)でマルチプルアライメントを行なって、結果を確認する。

mafft --auto input.fa > output.aln
alan output.aln

f:id:kazumaxneo:20180814121931j:plain

 

 引用

GitHub - mpdunne/alan: Alignment viewer for linux terminal

メタゲノムのシミュレータ CAMISIM

 16q rRNAアンプリコンとショットガンメタゲノムシーケンシングは、健康や病気に関するヒトマイクロバイオーム研究に広範に使われている[prepirntより ref.1, 2] 。私たちはその後、天然に存在する微生物群集は、生物多様性の広範な範囲をカバーしていることを学んだ。おそらくhalf a dozenのレベルから1万を超える微生物のpopulationsレベルの微生物多様性を含むことができ、代表的な分類群は大きく異なる可能性がある[ref.9-12]。これらの多様なコミュニティを分析することは困難である。この問題は、データ生成における幅広い実験設定の使用と、短期および長時間のシークエンシング技術の急速な進化によって悪化する[ref.13,14]。生成されるデータが非常に多様であるため、特定の実験設定に対する現実的なベンチマークデータセットを生成する可能性は、計算メタゲノミクスソフトウェアを評価するために不可欠である。

 メタゲノム解釈のクリティカルアセスメントのinitiative であるCAMIは、コンピュテーションメタゲノムソフトウェアの幅広く客観的なパフォーマンス概要を生成することを目的としたコミュニティの取り組みである[rwef.15]。 CAMIはベンチマーキングの課題を編成し、データの生成、ソフトウェアの適用、結果の解釈[ref.16]など、すべての面で標準の開発と再現性を促進する。
 ここでは、最初にCAMIの最初のチャレンジで使用されたシミュレートされたメタゲノムデータセットを生成するために書かれたCAMISIMについて説明する。 CAMISIMの有用性をいくつかのアプリケーションで実証する。著者らは、ヒトおよびマウスの腸内微生物のtaxonomyプロファイルから、複雑で多重反復のベンチマークデータセットを生成した[ref.1, 17]。数千もの小さな“minimally challenging metagenomes” をシミュレートし、ポピュラーなMEGAHIT [ref.18](紹介)やmetaSPAdes [ref.19]アセンブラ紹介)で、シーケンシングカバレッジ、ゲノムの進化的な相違、シーケンシングエラープロファイルなどの変化に伴う影響を特徴付けた。

 

 

f:id:kazumaxneo:20180810122531j:plain

CAMISIMのワークフロー。 Preprintより転載。

 

マニュアル

https://github.com/CAMI-challenge/CAMISIM/wiki/User-manual

 

CAMISIMに関するツイート。

インストール

依存

python 2.7.10

  • Biopython
  • BIOM
  • NumPy
  • Matplotlib

Genome annotation

  • Hmmer3 or RNAmmer 1.2(RNAmmer is a wrapper of Hmmer2. Hmmer uses hidden markov profiles to search marker genes in sequences.)
  • Mothur(A multi tool program. Alignment of sequences and clustering.)
  • MUMmer(A genome alignment software)  

Perl 5

Simulation

  • ART(ART is a set of simulation tools to generate synthetic next-generation sequencing reads.)
  • wgsim(Read simulator which offers error-free and uniform error rates.)
  • NanoSim(Read simulator for the generation of Oxford Nanopore Technologies (ONT) reads. )
  • PBsim(Read simulator for generating Pacific Biosciences (PacBio) reads.)
  • SAMtools 1.0

本体 Github

依存が多いのでここではdockerコンテナを使う。

docker pull cami/camisim:latest

> python metagenomesimulation.py -h

 docker run --rm -it -v /Users/user/Documents/docker_share/:/home/ cami/camisim metagenomesimulation.py -h

usage: python metagenomesimulation.py configuration_file_path

 

    #######################################

    #    MetagenomeSimulationPipeline     #

    #    Version 0.0.6                    #

    #######################################

 

    Pipeline for the simulation of a metagenome

 

optional arguments:

  -h, --help            show this help message and exit

  -v, --version         show program's version number and exit

  -silent, --silent     Hide unimportant Progress Messages.

  -debug, --debug_mode  more information, also temporary data will not be deleted

  -log LOGFILE, --logfile LOGFILE

                        output will also be written to this log file

 

optional config arguments:

  -seed SEED            seed for random number generators

  -s {0,1,2}, --phase {0,1,2}

                        available options: 0,1,2. Default: 0

                        0 -> Full run,

                        1 -> Only Comunity creation,

                        2 -> Only Readsimulator

  -id DATA_SET_ID, --data_set_id DATA_SET_ID

                        id of the dataset, part of prefix of read/contig sequence ids

  -p MAX_PROCESSORS, --max_processors MAX_PROCESSORS

                        number of available processors

 

required:

  config_file           path to the configuration file

ERROR: 0

——

 

> python genomeannotation.py

$ docker run --rm -it -v /Users/user/Documents/docker_share/:/home/ cami/camisim genomeannotation.py -h

usage: python genomeannotation.py configuration_file_path

 

#######################################

#    GenomeAnnotationPipeline         #

#    Version 0.0.6                    #

#######################################

 

Pipeline for the extraction of marker genes, clustering and taxonomic classification

 

optional arguments:

  -h, --help            show this help message and exit

  -v, --version         show program's version number and exit

  -verbose, --verbose   display more information!

  -debug, --debug_mode  tmp folders will not be deleted!

  -log LOGFILE, --logfile LOGFILE

                        pipeline output will written to this log file

 

optional config arguments:

  -p MAX_PROCESSORS, --max_processors MAX_PROCESSORS

                        number of available processors

  -s {0,1,2,3}, --phase {0,1,2,3}

                        

                        0 -> Full run (Default)

                        1 -> Marker gene extraction

                        2 -> Gene alignment and clustering

                        3 -> Annotation of Genomes

 

required:

  config_file           path to the configuration file of the pipeline

——

 

> python metagenome_from_profile.py

ヘルプなし。

 

ラン

ホストの指定パスをinputとoutputとして認識してランする。初回はNCBI database をダウンロードするので時間がかかる。

docker run -it -v "/Volumes/test/CAMISIM/defaults:/input:rw" \
-v "/Volumes/test/CAMISIM/out:/output:rw" \
cami/camisim metagenome_from_profile.py \
-p /input/mini.biom -o /output

 

 

 作成中。

 

引用

CAMISIM: Simulating metagenomes and microbial communities.

Adrian Fritz, Peter Hofmann, Stephan Majda, Eik Dahms, Johannes Droege, Jessika Fiedler, Till R. Lesker, Peter Belmann, Matthew Z. DeMaere, Aaron E. Darling, Alexander Sczyrba, Andreas Bremges, Alice C. McHardy

bioRxiv, 300970. doi:10.1101/300970 

 

参考

Critical Assessment of Metagenome Interpretation—a benchmark of metagenomics software

Alexander Sczyrb, Peter Hofmann, Alice C McHardy

Nature Methods volume 14, pages 1063–1071 (2017)

 

ターゲットキャプチャシーケンシングをシミュレートする capsim

 高スループットシークエンシング(HTS)は、費用対効果が高く時間効率の良いサンプルの完全な遺伝情報を得る能力を持ち、ゲノム研究に大きく革命をもたらした。多くの臨床応用において、作用可能な領域のパネルのみが調査対象である(Bellos et al、2014; Samorodnitsky et al、2015)(Diagnostic Panels.  e.g., Takara )。これらの分析では、研究者は、合成されたオリゴヌクレオチド(プローブ)のプールを用いて、ハイブリダイゼーションを用いて興味のあるゲノム断片を選択的に捕捉するターゲットキャプチャシーケンシングプロトコルを使用することが多い(Gnirke et al、2009)。これにより、ゲノムシーケンシング全体と比較して、より低いコストでより迅速な結果を得ることができ、臨床検査でのスケーラブルなアプローチが可能になる。

 計算シミュレーションは、HTSデータ分析ツールの開発とベンチマーキングに不可欠である(Escalona et al、2016)。In silicoでのシミュレーションデータは、実際のデータよりも安価に生産することができる。それらは制御された条件下で生成され、完全に特徴づけることができる。さらに、シミュレーションは、研究者がシーケンシングプロトコルの性能を評価し、実験を実施する前に設計を最適化するのに役立つ。whole genome sequencing(Escalona et al、2016)とtargeted exome sequencing (Kim et al、2013)には数多くのシミュレータが利用可能だが、私たち(著者ら)の知る限りでは、現在のところキャプチャプロセスのダイナミクスをin silicoでシミュレートできるツールは存在しない。このようなツールは、計算解析パイプラインだけでなく、キャプチャデザインの効率を評価するのにも役立つと考えている。

 この論文では、目標とするキャプチャシーケンシングデータをシミュレートする必要性を満たすためのソフトウェアパッケージであるCapSimを紹介する。プローブセットが与えられると、CapSimは、in silicoにおけるプローブハイブリダイゼーションのダイナミクスをシミュレートし、シーケンシングされるフラグメントセットを生成する。既存のHTSシミュレーションツールとは異なり、CapSimは、フラグメンテーション、フラグメントキャプチャ、シーケンシングなど、シーケンス処理のあらゆる段階をエミュレートする。 CapSimはJavaで書かれており、あらゆるコンピューティングプラットフォーム上でネイティブに実行できるため、バイオインフォマティクスコミュニティに簡単にアクセスできる(どのような環境のバイオインフォマティシャンでも利用できる)。

 シーケンシングプロセスはDNAの断片化から始まる。CapSimは所定のサイズ分布を有するゲノム配列からフラグメントを繰り返しサンプリングしてこのプロセスをシミュレートする。いくつかのデータセットからの断片サイズ分布に最もよく適合することがわかったlog-logistic分布を用いて断片長をモデル化する(論文の補足資料参照)。

 ターゲットキャプチャシーケンシングの次のステップでは、DNA断片にハイブリダイゼーションにした標的プローブとそのDNA断片はビーズからpull downされる。ビーズは標的プローブとハイブリダイズしたDNA断片と特異的に結合している。CapSimはプローブをDNA断片にマッピングすることでこのプロセスをシミュレートする。計算上効率的であるために、CapSimは最初にプローブをゲノム配列にマップし、DNA断片がゲノムからサンプリングされると、greedy algorithmを使用してDNA断片に結合できるプローブの最大数を決定する。我々(著者ら)は、キャプチャされる確率がバインドしたプローブの数に比例し、断片長に反比例するキャプチャプロセスの確率論的性質をモデル化した。ハイブリダイゼーションのダイナミクスのこのシミュレーションは、Kim et al. (2013)よりもより現実的なキャプチャシーケンシングデータを生成することが示されている。 Kim et al. (2013)は、キャプチャシークエンシングをシミュレーションするための唯一のツールである(論文の結果参照)。

 次いで、キャプチャされたフラグメントはin silico sequencingされる。 Illuminaシーケンシングのために、CapSimはクラスター形成のDNA断片分布を導入ししている。 それはLog-Logistic分布を使用して、キャプチャステップからシミュレートされたDNA断片からサンプリングを行う(論文の補足資料を参照)。これらの選択されたDNA断片は、DNA断片の両末端からエラーありで配列が読み取られペアエンドシークエンシングシミュレーションに使用される。 PacBioシークエンシングのために、CapSimは、所与の分布からのポリメラーゼが読む長さをシミュレートする(補足資料参照)。DNA断片トのシークエンシングをシミュレートする際、CapSimはポリメラーゼが読める長さに達するまで、PacBioエラープロファイルを使い、2つの鎖のコピーを交互に行う。

 

 

BioinfomaticsのSEQUENCE ANALYSISに発表された論文です。大半のデータはsupplementaryにあります。

 

capsimに関するツイート。

 

Japsa(Just Another JAva Package for Sequence Analysis)マニュアル

http://japsa.readthedocs.io/en/latest/install.html

 

ダウンロード

mac os 10.13、java1.6環境でテストした。

Github

リリースよりjapsaをダウンロードする。

https://github.com/mdcao/japsa/releases

tar zxvf JapsaRelease.tar.gz
cd JapsaRelease
./install.sh

コンソールに表示される指示に従いインストールする。最後の質問のHDF libraryがなければ"none"と打つ。ここではbin/にパスを通した。

cd bin/
echo 'export PATH=$PATH:`pwd`' >> ~/.bash_profile && source ~/.bash_profile

jsa.sim.capsim --help

$ jsa.sim.capsim --help

Simulate capture sequencing

 

Usage: jsa.sim.capsim [options]

Options:

  --reference=s   Name of genome to be 

                  (REQUIRED)

  --probe=s       File containing probes mapped to the reference in bam format

                  (default='null')

  --logFile=s     Log file

                  (default='-')

  --ID=s          A unique ID for the data set

                  (default='')

  --miseq=s       Name of read file if miseq is simulated

                  (default='null')

  --pacbio=s      Name of read file if pacbio is simulated

                  (default='null')

  --fmedian=i     Median of fragment size at shearing

                  (default='2000')

  --fshape=d      Shape parameter of the fragment size distribution

                  (default='6.0')

  --num=i         Number of fragments 

                  (default='1000000')

  --pblen=i       PacBio: Average (polymerase) read length

                  (default='30000')

  --illen=i       Illumina: read length

                  (default='300')

  --seed=i        Random seed, 0 for a random seed

                  (default='0')

  --help          Display this usage and exit

                  (default='false')

 

——

 

 

ラン 

genepanelのfasta配列が見つからなかったので、biomartで生成する。ここでは、指定した領域から配列を出力している(gene idを持っているなら、biomartでidから配列を生成する方が王道。IDを別種のID、または遺伝子名に変換したいならDavidを使う)。

 

1、配列の準備

Gene panelの領域は、Perkin Elmer: Amplicon Panel Kitsの領域を使う。

http://www.biooscientific.com/Illumina-Amplicon-Panels

ここではColorectal Cancer -1 kitのターゲット領域を選んだ。 

http://www.biooscientific.com/NEXTflex-Colorectal-Cancer-1-Amplicon-Panel

 

Ensembl - BIomartにアクセスする。

http://asia.ensembl.org/biomart/martview/f5938f5b76d2d943d8c4a176709e376a

f:id:kazumaxneo:20180810003230p:plain

Ensembl Genes 93を選択、Human genome GRCh38を選択。

 

左のタブからFilterを選択。様々な条件で配列などを抽出できる。代表的なデータベースのIDからも配列は抽出できるが、ここではbedファイルしかないので、ポジション情報から配列を抽出する。Multiple regionにbedの配列をペーストする。bedファイルのタブはあらかじめ":"に置換してある。

f:id:kazumaxneo:20180810003114p:plain

 貼ったところ。

 

左のメニューからAttributeタブを選択。SequencesとcDNA sequencesを選択。

f:id:kazumaxneo:20180810003900p:plain

SequencesとcDNA sequencesを選んだところ。

 

左上のResultsをクリックして、配列を生成する。Compressed fileをダウンロードする。

f:id:kazumaxneo:20180810004054p:plain

ターゲット領域のMulti-FASTAファイルが準備できた。テストランでは、これをprobes.faとして使う。 capsimのランに移る。

 

2、capsim実行。

probes.faをリファレンスにmappingする。GRCh38はEnsembl humanからダウンロードした。

bowtie2-build ref.fa ref
bowtie2 --local --very-sensitive-local --mp 8 --rdg 10,8 --rfg 10,8 -k 10000 -f -x ref -U probes.fa -S probes.sam
samtools sort -O BAM probes.sam > probes.bam && samtools index probes.bam

capsimを実行。

#pacbio
jsa.sim.capsim --reference ref.fa --probe probes.bam --ID someid --fmedian 5000 --pacbio output --pblen 3000 --num 20000000

#Miseq
jsa.sim.capsim --reference ref.fa --probe probes.bam --ID someid --fmedian 500 --miseq output --illen 300 --num 20000000
  •  --ID    A unique ID for the data set  (default='')
  • --fmedian    Median of fragment size at shearing (default='2000')
  • --pacbio    Name of read file if pacbio is simulated (default='null')
  • --pblen      PacBio: Average (polymerase) read length   (default='30000')
  • --num        Number of fragments  (default='1000000')
  • --miseq     Name of read file if miseq is simulated (default='null')
  • --illen         Illumina: read length  (default='300')

マッピング結果。

figv -g ref.fa miseq.bam,Colorectal_Cancer_-1.bed

f:id:kazumaxneo:20180810090305p:plain

 

引用

Simulating the dynamics of targeted capture sequencing with CapSim
Cao MD, Ganesamoorthy D, Zhou C, Coin LJM

Bioinformatics. 2018 Mar 1;34(5):873-874. 

 

参考

TruSeq Amplicon Cancer Panel

https://jp.illumina.com/content/dam/illumina-marketing/apac/japan/documents/pdf/2012_illumina_truseq_amplicon_cancer_panel.pdf

 

 

 

 

 

YSTRなどのショートタンデムリピートを探す STRScan

 

 マイクロサテライトまたは単純配列反復(SSR)とも呼ばれる短いタンデムリピート(STR)は、タンデム反復ユニット(1〜6 bps)を約2〜30個含む短いストレッチのDNAである。 STRは、ヒトなどの哺乳動物ゲノムを含む多くの原核生物および真核生物ゲノムに存在する[論文より 1,2]。 50万以上のSTRがヒトゲノムで特徴づけられており、ヒトゲノム全体の約3%を構成している[ref.3]。多型性(polymorphism)が高いため、STRは遺伝マーカーとして一般的に使用されている[ref.4-7]。特に、STR遺伝子座の小さなセットは、1つの(時には未知の)供給源からの少量のヒトDNAにおいてPCRを用いて複数のSTR遺伝子座が増幅された同一性および親の検査に使用することができる(ref.8,9)。 PCR産物の長さは、他の供給源由来の1つ以上のヒトDNAサンプル(例えば、フォレンジックデータベース中)と比較される。このSTRタイピング手続きは、大部分が標準化されており、そのようなテストの対象となるSTR遺物は、STRBase [ref.10]などのpublucデータベースに収集されている。

 STRは大部分が「ジャンクDNA」と考えられているが、STRの中にはタンパク質コード遺伝子が含まれており、その産物は転写調節に関与するグルタミンリッチなドメインなど高等生物に機能的役割を果たすことが示されている[ref.11]。非コード領域のSTRでも、それらの下流遺伝子の発現調節に関与している可能性がある[ref.12]。特に、トリヌクレオチドリピートの不安定な拡大は、ヒト疾患に関連することが知られている[ref.13]。抜群の例はハンチントン病であり、ハンチントン(HTT)遺伝子におけるCAGトリプレットのタンデムリピートの拡大に​​よって引き起こされ、脳変性に至る可能性のある異なるタンパク質形態をもたらす遺伝的な神経変性疾患である[ref.14]。したがって、疾患感受性対立遺伝子におけるSTRプロファイリングは、これらの遺伝的障害を遺伝する高いリスクを有する個体のための遺伝子検査ツールとしてしばしば用いられる[ref.15]。

 STRの伝統的な実験分析は、プライマーをSTRに隣接するユニークな領域に設計し、PCRによる標的STR遺伝子座を増幅し、続いてゲル電気泳動を用いたPCR産物の長さ測定が行われる。長さはSTRのコピー数によって決まる。近年、次世代シーケンシング(NGS)技術の急速な進歩により、全ゲノムシーケンシング(WGS)がより手頃なものになっている。Tandem repeat finder[ref.16]のような従来のソフトウェアツールは、ヒトゲノムのようなアセンブリされたゲノム配列から新規なSTRを検出することができる[ref.17]。 WobデータのSTRプロファイリングに直接適用できるlobSTR [ref.18]やSTR-FM [ref.19]などの新しいソフトウェアツールとパイプラインも開発されている。最近の研究では、NGSデータからSTRを分析する能力が示されており、ヒト染色体STR(Y-STR)のプロファイルの分析を通じて、ヒトゲノム配列データからヒト個体の姓を推測することができ、オンラインの系譜(genealogy)データベース[ref.20]などがある。ゲノムワイドのSTRプロファイリングツールは、集団におけるSTR変異の調査を可能にしている[ref.19,21,22]。かなりの数のSTR遺伝子座がヒト集団中に広範に発現され、これはヒトゲノムにおける新規なregulatory variantsセットを表すかもしれないことが示された(ref.23)。

 本稿では、次世代シーケンシング(NGS)データのSTRプロファイリングのためのスタンドアロンのソフトウェアツールSTRScanを紹介する。ここでは、STRプロファイリングのためにターゲットアプローチを採用した。ゲノム全体スケール(lobSTRまたはSTR-FMの目標)ですべてのSTRを採掘するのではなく、ユーザーが定義したSTR locusのサブセット、つまり法医学または遺伝子検査に有用であり[24]、よって時間がかかるゲノムワイドのマッピング手順を回避するということである。結果として、本発明者らの方法は、リファレンスゲノム中から線状DNA配列として表されるSTR遺伝子座を配列比較することによる制限は起きず、DNA配列中のSTR同定のためファインチューニングしたアラインメントアルゴリズムを採用できることになる。本発明者らの方法は、全ゲノムシーケンシングデータからのマイニングに加え、PCRによるSTRエンリッチメントを行ったターゲットSTRのNGSデータに直接適用することができる。

 STRScanにおいて、各STR遺伝子座は1つまたは複数のタンデムコピーを含むパターンとその間にある配列によって表される。これはリファレンスゲノムから構築できる(例えばヒト)、そしてgreedy seed-extension戦略を使用して、シーケンシング配列の各STR遺伝子座を同定する(一部略)。

 著者らは、Sangerシーケンサー[ref.26]およびIlluminaシーケンサー(1000 Genomes Project [ref.27]によって生成された)の全ゲノムシーケンシング(WGS)データに対してSTRScanを試験した。 STRScanは、lobSTRやSTR-FMのような既存のソフトウェアツールと比較して、同等かそれ以下の計算時間を使用しながら、NGSデータからより多くのSTRを大幅に(平均20%)識別することができた。STRScanはDNA増幅と続く次世代シーケンシングによるSTRタイピングに使う準備が整っている。

 

STRScanに関するツイート。

 

STRScan HP

f:id:kazumaxneo:20180807213227p:plain

 

インストール

cent OS6のpython2.7.10でテストした。

依存

  • python2

プログラム本体は公式HPからダウンロードする。

http://darwin.informatics.indiana.edu/str/

tar -zxvf STRScan.tar.gz
cd STRScan/
make #binaryがはあるがsegmentation errorを起こすなら再ビルド

./STRScan

$ ./STRScan -h

./STRScan: option requires an argument -- 'h'

STRScan -i inpfile -p patfile -o outfile -g ngap -w window

-i InpFile: The input file name of sequence reads

-p PatFile: The input file name of short tandem repeat (STR) patterns

-o OutFile: The output file with profiled repeat copies 

-g maximum allowed gaps

-w maximum allowed windows of copy numbers

-q: input file is in fastq format (default in fasta format) 

-l: minimum lengths of the supported alignments of the spanning sequence

 

 

テストラン

cd STRScan/example/

 ランにはshort tandem repeatファイルが必要になる。テストではYSTR_pat(Y染色体のSTR:wiki)を使う。

> head -n 20 YSTR_pat

$ head -n 20 YSTR_pat

>DYS19 (TCTA)12CCTA(TCTA)3

atgccacccttttattatttctacggatattacttggactggaagacaag

gactcaggaatttgctggtcaatctctgcacctggaaatagtggctgggg

caccaggagtaatacttcgggccatggccatgtagtgaggacaaggagtc

catctgggttaaggagagtgtcactata

---

aaacactatatatatataacactatatatataatactatatatatattaa

aaaacactataacagaaactcagtagtcatagtgaaatcaaaaaataatc

acagtcaatttgatctcatacctagactgaaatatgaaacttcaaaagaa

aagaatgttaagaactttgggcttgtcaaaattttcctacatagataa

---

>DYS19  (TAGA)3TAGG(TAGA)12

gattattttttgatttcactatgactactgagtttctgttatagtgtttt

ttaatatatatatagtattatatatatagtgttatatatatatagtgttt

---

tatagtgacactctccttaacccagatggactccttgtcctcactacatg

gccatggcccgaagtattactcctggtgccccagccactatttccaggtg

---

>DYS385a (GAAA)14

agtgcatgtaatcccagctacttgggaggctgaggcagggtaattgtttg

このshort tandem repeatリストが指定したFASTAファイルtest.faに見つかるかを調べる。

STRScanを実行する。

STRScan -i test.fa -p YSTR_pat -o test-5.out -g 5 -w 3 -l 5

検索結果のFASTAとhit部位の情報ファイルが出力される。

 

残り3つのshort tandem repeatファイルでも同様にランする。

STRScan -i test.fq -p YSTR_pat -o test_fq-5.out -g 5 -w 3 -l 5 -q
STRScan -i test.fq -p CODIS_pat -o test_c-5.out -g 5 -w 3 -l 5 -q
STRScan -i test.fq -p CODIS_pat -o test_c_fq-5.out -g 5 -w 3 -l 5 -q

bedに変換する。

python ../ParseResult.py -r1 ref_CODIS.bed -f1 test_c-5.out | awk '!NF || !seen[$0]++' > test_c-5.bed

> cat test_c-5.bed

$ cat test_c-5.bed

CSF1PO -8|1 4 11 13

D7S820 NA|0 4 0 13

  1. First column: STR Marker
  2. Second column: (STRScan allele -Reference allele)*STR period | Number of supporting reads
  3. Third column: STR period Fourth column: STRScan reported allele
  4. Fifth column: reference allele  

 

ラン 

fastqの分析。fastaフォーマットの場合は"-q"を外す。

STRScan -i input.fa -p pattern_file -o output -g 5 -w 3 -l 5 -q
  • -i    the reads file (in fasta or fastq format; default is fasta, for fastq input, use -q)
  • -p   the input file of the STR library (in the format as you know)
  • -o   there output report file (there is also another file with .seq ext with the matched reads sequences in fasta format)
  • -g   alignment bandwidth, default 5
  • -l    the minimum length for the spanning sequence (upstream and downstream), default 5
  • -w   the maximum copy number exceeding the given copy number of STR, default 3 (e.g., for (CTT)10, the algorithm will allow for maximum 13 copies of the repeat units).
  • -q    input file is in fastq format (default in fasta format) 

出力はテストと同じようなbedファイルと、サポートするリード情報になる。

 

引用

STRScan: targeted profiling of short tandem repeats in whole-genome sequencing data

Haixu Tang and Etienne Nzabarushimana

BMC Bioinformatics. 2017; 18(Suppl 11): 398.