macでインフォマティクス

macでインフォマティクス

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

アンプリコンシーケンスのバリアントコーラー UNDR ROVER

 

 超並列シーケンシング(MPS)[論文より ref.1]のためにHi-Plex(www.hiplex.org)と呼ばれる高度に多重化されたPCRベースのターゲット濃縮システムを開発した。 Hi-Plexは、非常に正確な結果を得ることができるシンプルで低コストのプロトコルである。その主要な特徴の1つは、サイズ選択によるオフターゲット増幅の除去を容易にする均一なライブラリーサイズを定義する能力であり、ペアエンドシーケンシングと組み合わせて、各アンプリコンについてのペアエンドの完全なオーバーラップを可能にする。後者は、バリアントの検出およびシーケンスエラーによって引き起こされるアーティファクトのフィルタリング両方において高度な厳密さを可能にする。以前、Hi-Plexシーケンス [ref.2](参考ページ)によって生成されたオーバーラップリードを使うariable CallingツールであるROVERを開発し、PALB2のコード領域の遺伝的な変異のスクリーニングにHi-PlexとROVERを適用した。結果は、以前のスクリーニングによって同定された全60変異を全て検出し、偽陽性はゼロだった[ref.3、4]。

 ROVERは、ターゲットとしてアンプリコン領域のゲノム座標をタブ区切り形式で記述したファイルと、リファレンスゲノムにマップされたペアエンドリードを含むBAM形式[ref.5]の1つ以上の配列ファイルを必要とする。出力はVCF形式のバリエーションのリストである[ref.6]。 ROVERは、SNVおよびsmall indelを検出することができる。バリアントは、オーバーラップするペアエンドリードの両方で同じ位置に表示される場合にのみ、ROVERによってレポートされる。

 ROVERを用いてバリアントを検出するのに最も計算量が高価になる部分は、リファレンスゲノムへのマッピング(またはアライメント)である。以下で説明する1つの指標実験では、Bowtie(http://bowtie-bio.sourceforge.net/)[ref.7]とのマッピングに要した時間は、ROVERバリアントコールパイプライン全体の約78%を構成していた。リードのマッピングは全エキソームおよび全ゲノムDNAシークエンシングパイプラインの標準的な部分であるが、これまでにAmplivar [ref.8]によって実証されているように、これはmassively parallel sequencing (MPS) ベースのPCRによって避けることができる。これは、各シーケンスリードの5 '末端がPCR増幅のプラ​​イマー配列で始まり、そして既知の配列を増幅するからである。後者の情報は、プライマー設計中に決定される。プライマーペアが介在する配列のゲノム座標を同定するのでリードはリファレンスにマップする必要はない。オプションで、ペアの少なくとも1つのリードのプライマーに続く小さな配列がリファレンスと同一であることを確認することで、正しい位置にマッピングされるという確信をさらに高めることもできる。

 特定のリードの開始座標を決定したら、そのシーケンスをリファレンスと比較することができる。挿入と削除検出を可能にするためには、Needleman-Wunschアルゴリズム[ref.9]のスタイルでギャップアライメントを実行する必要がある。しかしながら、このアルゴリズムの複雑さは、アライメントされたシーケンスの長さの2次関数であり、したがって、入力における各リードについて計算するコストが大きくなる。幸運なことに、入力のほとんどのリードはリファレンスと同じになるか、少数の不一致でしか違いがないため、このコストを避けることができる。このような場合、リードとリファレンスの単純な線形比較で十分である。線形の比較が失敗した場合にのみ、ギャップアライメントアルゴリズムに戻る。

 Amplivarは、UNDR ROVERの基礎となる構想に似ているが、異なる仕組みを適用している。 Amplivarはプライマー配列を使用してアンプリコンとリードを関連させ、リードをグループとしてアライメントさせることで計算上のオーバーヘッドを削減する(BLAT [ref.10]を使用)。さらに、Amplivarは(SeqPrep [ref.11]を使用して)オーバーラップリードをマージするが、UNDR ROVERは両方のリードを保持して厳密なフィルタリングシステムの一部として一致するかテストする。

 UNDR ROVERは、Hi-Plexターゲットシーケンスをサポートするように設計されているが、シーケンシングリードで遺伝子特異的プライマー配列を保持する他のアンプリコンベースのターゲットシーケンスシステムとも互換性がある。 AmpliSeqのデータは、例えば、遺伝子特異的プライマーがライブラリー生成の間に大部分切断されるため、適合性がない。重要な考慮点は、UNDR-ROVERは、ペアエンドリードのかなりのオーバーラップを示すシーケンシングデータを処理することを意図しているため、これを達成していないシステムでの使用はお勧めしない。

 

インストール

mac os 10.13のvirtualenv環境でテストした。

Github

 公式ではvirtualenvを使い、仮想環境で導入、実行することを勧めている。virtualenvを使ったインストールは以下のように行う。

virtualenv --system-site-packages undr_rover_dev source
undr_rover_dev/bin/activate
pip install git+https://github.com/bjpop/undr_rover

> undr_rover --version

$ undr_rover --version

undr_rover 2.0.0

undr_rover

$ undr_rover

usage: undr_rover [-h] [--version] --primer_coords PRIMER_COORDS

                  --primer_sequences FILE

                  [--primer_prefix_size PRIMER_PREFIX_SIZE]

                  [--kmer_length KMER_LENGTH]

                  [--kmer_threshold KMER_THRESHOLD]

                  [--primer_bases PRIMER_BASES] [--proportionthresh N]

                  [--absthresh N] [--qualthresh N] [--overlap OVERLAP]

                  [--max_variants N] --reference FILE [--id_info FILE] --out

                  FILE [--log FILE] [--coverdir COVERDIR] [--fast]

                  [--snvthresh N] [--genotype] [--ploidy {1,2}]

                  [--error ERROR]

                  fastqs [fastqs ...]

undr_rover: error: too few arguments

(undr_rover_dev) kamisakBookpuro:undr_rover_dev kamisakakazuma$ undr_rover -h

usage: undr_rover [-h] [--version] --primer_coords PRIMER_COORDS

                  --primer_sequences FILE

                  [--primer_prefix_size PRIMER_PREFIX_SIZE]

                  [--kmer_length KMER_LENGTH]

                  [--kmer_threshold KMER_THRESHOLD]

                  [--primer_bases PRIMER_BASES] [--proportionthresh N]

                  [--absthresh N] [--qualthresh N] [--overlap OVERLAP]

                  [--max_variants N] --reference FILE [--id_info FILE] --out

                  FILE [--log FILE] [--coverdir COVERDIR] [--fast]

                  [--snvthresh N] [--genotype] [--ploidy {1,2}]

                  [--error ERROR]

                  fastqs [fastqs ...]

 

Find variants from fastqs via a mapping-free approach.

 

positional arguments:

  fastqs                Fastq files containing reads.

 

optional arguments:

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

  --version             show program's version number and exit

  --primer_coords PRIMER_COORDS

                        Primer coordinates in TSV format.

  --primer_sequences FILE

                        Primer base sequences as determined by a primer

                        generating program.

  --primer_prefix_size PRIMER_PREFIX_SIZE

                        Size of primer key for blocks in dictionary.

  --kmer_length KMER_LENGTH

                        Length of k-mer to check after the primer sequence.

                        Defaults to 30.

  --kmer_threshold KMER_THRESHOLD

                        Number of single nucleotide variants deemed acceptable

                        in kmer.

  --primer_bases PRIMER_BASES

                        Number of bases from primer region to use in gapped

                        alignment.Helps with variant calling near the edges of

                        a block. Defaults to 5.

  --proportionthresh N  Keep variants which appear in this proportion of the

                        read pairs. For a given target region, and bin

                        otherwise. Defaults to 0.05.

  --absthresh N         Only keep variants which appear in at least this many

                        read pairs. Defaults to 2.

  --qualthresh N        Minimum base quality score (phred).

  --overlap OVERLAP     Minimum proportion of block which must be overlapped

                        by a read. Defaults to 0.9.

  --max_variants N      Ignore reads with greater than this many variants

                        observed. Defaults to 20.

  --reference FILE      Reference sequences in Fasta format.

  --id_info FILE        File containing rs ID information in VCF format.

  --out FILE            Name of output file containing called variants.

  --log FILE            Logs progress in specified file, defaults to stdout.

  --coverdir COVERDIR   Directory to write coverage files, defaults to current

                        working directory.

  --fast                Use gapped alignment less often, leading to faster run

                        time.

  --snvthresh N         Distance between two single nucleotide variants before

                        going to a gapped alignment.

  --genotype            Compute genotypes for SNVs. Defaults to False.

  --ploidy {1,2}        Ploidy for genotyping 1 = haploid, 2 = diploid.

                        Defaults to 2.

  --error ERROR         Expected base read error rate. Defaults to 0.002.

解析が終わったら"deactivate"で仮想環境を終了する。 

 

ラン

テストラン。

git clone https://github.com/bjpop/undr_rover 
cd undr_rover/example/

undr_rover --primer_coords example_coords.txt --primer_sequences \ example_primers.txt --reference hg19.fa --out example.vcf \
--genotype example_R1.fastq example_R2.fastq

出力はVCF形式に準じている。

 

 

引用

UNDR ROVER - a fast and accurate variant caller for targeted DNA sequencing.

Park DJ, Li R, Lau E, Georgeson P, Nguyen-Dumont T, Pope BJ

BMC Bioinformatics. 2016 Apr 16;17:165.