次世代シークエンシング（NGS）技術は、以前はあまりにも労働集約的であると考えられていた実験を日常的な作業に変えた（Metzker、2010）。 NGSの1つの用途は、genetic diversityを定量化するために遺伝的にheterogousなpopulationsのシーケンシングである。Genetic diversityは、例えば、HIVおよびHCVなどのRNAウイルスによる感染において主な関心事である。これらのシステムでは、病原体の高い突然変異率と病原体間の組換えとの組み合わせにより、ウイルス性疑似種（quasispecies）と呼ばれる異なるが関連する個体集団が生じる。異なるウイルス粒子が単一の細胞に感染すると、後の事象として組換えが起こり得る。我々（著者ら）は、この集団における異なるウイルス株をハプロタイプとして示す。ウイルス準種の特徴を研究することで、宿主における病原体の進化のメカニズムを明らかにすることができる。例えば、準種の多様性は病原性（Vignuzzi et al、2006）、immune escape （Nowak et al、1991）および薬剤耐性（Johnson et al、2008）に影響を及ぼすことが示されている。
The quasispecies equationは、突然変異選択プロセスに従って進化するRNAウイルス集団の数学モデルである（Eigen、1971）。モデルのダイナミクスは、複製時にあるウイルスのハプロタイプ（または株）を別のものへと変換する突然変異と、異なる株の複製速度の変動を説明する選択用語とによって記述される。 HIVおよびHCVを含む多くの臨床的に関連性の高いウイルスでは組換えが頻繁であることが知られているが、突然変異プロセスは一般に点突然変異の結果としてのみ考慮される。例えば、HIVの組換え率は、その点突然変異率より約10倍高いと推定される。したがって、quasispecies modelは、突然変異および組換えの両方を説明するために拡張されている（Boerlijst et al、1996）。平衡状態では、モデルは、ウイルス集団が1つまたは少数のハプロタイプによってdominatであり、、常に生成される低頻度の突然変異株の雲に囲まれていると予測する。
最近のNGS技術は、1回の実験で数百万のDNAシーケンシングリードを生成することによって、前例のないレベルでウイルスquasispeciesを観察することを可能にする。しかしながら、この高い収率はコストがかかる。リードは、通常、最新の技術では700bpまでの短差で、最小のウイルスゲノムよりずっと小さく、サンプル調製やシーケンシングエラーが原因でエラーが発生しやすい（Gilles et al、2011）。結果として、得られたデータは不完全でノイズが多いので、NGSを用いてウイルス集団を有意に特徴づけるには、シーケンシングデータを注意深く解析する必要がある（Beerenwinkel and Zagordi、2011）。
この論文では、突然変異と組換えを明示的にモデル化することにより、NGSデータに基づいてウイルスquasispeciesを推定することを目指している。我々（著者ら）は、ウィルス集団、すなわち、ハプロタイプ分布およびNGSによるプロービングを生成するために、隠れマルコフモデル（HMM）を使用する。本発明者らのモデルでは、ハプロタイプは、ハプロタイプがどの配列から派生するかを選択する、HMMの状態の切り替えとして記述される少数の生成配列、および生成配列の位置特異的確率表によって記載される突然変異から生じる。シークエンシングリードは、観察エラーの対象となるハプロタイプから得られる。
本論文では、ウイルスquasispeciesの推定を行うための、すなわち患者内のウイルスのハプロタイプの分布を推定するための新規の確率的モデルを提示する。具体的には、実際の遺伝的多様性は、突然変異および組換えにより、generatorsと呼ばれる、いくつかの配列によって生成され、観察された多様性は追加のシーケンシングエラーに起因すると仮定する。我々（著者ら）は局所ハプロタイプ推論のモデルを提示する。すなわち、個々のリードによってカバーすることができる大きさのゲノム領域の集団構造を推定することを目的とするが、このモデルをグローバルなハプロタイプ推論、すなわちより長いゲノム領域に拡張することを目的とする。一般に、ローカル推論は、多くのアプリケーションにとって、より信頼性が高く、十分である。例えば、抗レトロウイルス治療の重要な標的であるHIVプロテアーゼ遺伝子は297bpであり、臨床診断のための共通のパイロシーケンシングプラットフォームであるRoche / 454 GS Juniorシーケンサーで400bp以上のリード値を得ることが現在の標準である（論文執筆時点）。局所的なハプロタイプの再構成はまた、全体的な再構築の出発点として使用することもできる（Astrovskaya、2011; Eriksson et al、2008; Prabhakaran et al、2010; Prosperi et al、2011; Zagordi et al、2011）。本モデルは、ウイルス性ハプロタイプの分布を高い信頼性で推定することができることを、ground truthにアクセスできるシミュレートされたデータに適用することで示し、そしてHIV感染患者から得たリードへの適用を示す。
QuasiRecombに関するツイート
特徴（Githubより）
- First algorithm that models the recombination process
- Allows position-wise mutation events
- Infers a parametric probability distribution from the underlying viral population
- Error correction by estimating position-wise sequencing error-rates
- Local, gene- and genome-wide reconstruction
- Reports SNV (single nucleotide variant) posteriors
- Incorporates paired-end information
- Uses PHRED scores to weight each base of each read
- Input may contain paired-end and single reads
- Supports reads of all current sequencing technologies (454/Roche, Illumina and PacBio)
- Suitable for amplicon and shotgun sequencing projects
- Reports reconstructed haplotypes and their relative frequencies
- Reports translated proteins in all three reading frames with their relative frequencies
- Input data can be in BAM or SAM format
インストール
本体 Github
リリースよりビルドされたjavaプログラムをダウンロードする。
> java -jar /Users/user/Downloads/QuasiRecomb.jar
$ java -jar /Users/user/Downloads/QuasiRecomb.jar
No input given
QuasiRecomb version: 1.2
Get latest version from http://bit.ly/QuasiRecomb
USAGE: java -jar QuasiRecomb.jar options...
-------------------------
=== GENERAL options ===
-i INPUT : Alignment file in BAM or SAM format.
-o PATH : Path to the output directory (default: current directory).
-K INT or INT-INT : The interval or fixed number of sequence generators, i.e. 1-4 or 2
In a grid enviroment the $SGE_TASK_ID.
In case of no input, K will be incremented as long as max BIC has not been reached, but will stop at K=5.
-m INT : The number of EM restarts during model selection (default: 5).
-t INT : The number of EM restarts for best K to find optimum (default: 50).
-r INT-INT : Only reconstruct a specific region.
-noRecomb : Do not allow recombination.
-quality : Account phred quality scores (slower runtime).
-printAlignment : Save alignment.txt in a human readable format.
-sampleProteins : Sample full-length protein sequences in three reading frames.
-coverage : If your dataset only contains a single region of interest,
regions with a minimum coverage of 100x, 500x, 1,000x and 10,000x are reported.
-refine : Can only be used after QuasiRecomb has been executed once before on the same dataset in the same directory.
Thins the number of haplotypes.
-noGaps : Ignore gaps; useful if data is 454 and gaps are only technical errors.
-conservative : Use this if the major haplotypes are only of interest.
-maxDel INT : Remove reads with more than INT consecutive deletions.
-maxPercDel DOUBLE : Remove reads with more than DOUBLE ratio of deletions, between 0.0 - 1.0
-unpaired : If read names are not unique and reads are single-end, prevent pairing and merging.
-------------------------
=== Technical options ===
-XX:NewRatio=9 : Reduces the memory consumption (RECOMMENDED to use).
-Xms2G -Xmx10G : Increase heap space.
-XX:+UseParallelGC : Enhances performance on multicore systems.
-XX:+UseNUMA : Enhances performance on multi-CPU systems.
-------------------------
=== EXAMPLES ===
java -XX:NewRatio=9 -jar QuasiRecomb.jar -i alignment.bam
java -XX:NewRatio=9 -jar QuasiRecomb.jar -i alignment.bam -conservative
java -XX:NewRatio=9 -jar QuasiRecomb.jar -i alignment.bam -K 1:10
java -XX:NewRatio=9 -jar QuasiRecomb.jar -i alignment.bam -noRecomb -r 790-2292
java -XX:+UseParallelGC -Xms2g -Xmx10g -XX:+UseNUMA -XX:NewRatio=9 -jar QuasiRecomb.jar -i alignment.bam
-------------------------
For further information, see http://bit.ly/QuasiRecomb-howto
-------------------------
実行方法
java -jar QuasiRecomb.jar -i alignment.bam
テスト
論文と同じSRR069887のデータを使いテストした。リファレンスはHXB2（リンク）を使用。
マッピング結果
ターゲット領域のリードのみ取り出す。
samtools view -@ 10 -b sort.bam chr:6321-6499 > target.bam
QuasiRecombのラン
java -jar QuasiRecomb.jar -i target.bam
> head -n 50 support/K1-snvs.txt
head -n 50 support/K1-snvs.txt
A C G T -
0 0 0 0 1 0
1 1 0 0 0 0
2 0 0 1 0 0
3 0 0 0 1 0
4 1 0 0 0 0
5 1 0 0 0 0
6 1 0 0 0 0
7 1 0 0 0 0
8 1 0 0 0 0
9 0 0 0 1 0
10 0 0 0 1 0
11 0 0 1 0 0
12 1 0 0 0 0
13 1 0 0 0 0
14 0 1 0 0 0
15 1 0 0 0 0
16 0 0 0 1 0
17 1 0 0 0 0
18 0 0 1 0 0
19 0 0 1 0 0
20 0 0 0 0 1
21 0 0 0 0 1
22 1 0 0 0 0
23 1 0 0 0 0
24 0 0 0 1 0
25 1 0 0 0 0
26 0 0 1 0 0
27 0 1 0 0 0
28 1 0 0 0 0
29 0 1 0 0 0
30 0 1 0 0 0
31 1 0 0 0 0
32 0 1 0 0 0
33 0 1 0 0 0
34 1 0 0 0 0
35 1 0 0 0 0
36 0 0 1 0 0
37 0 0 1 0 0
38 0 1 0 0 0
39 1 0 0 0 0
40 1 0 0 0 0
41 1 0 0 0 0
42 0 0 1 0 0
43 1 0 0 0 0
44 0 0 1 0 0
45 1 0 0 0 0
46 0 1 0 0 0
47 1 0 0 0 0
48 0.23398 0 0.76601 0 0
K1-summary.html
quasispecies.fasta
サマリープロット
R CMD BATCH support/coverage.R
R CMD BATCH support/modelselection.R
引用
Probabilistic Inference of Viral Quasispecies Subject to Recombination
Armin Töpfer, Osvaldo Zagordi, Sandhya Prabhakaran, Volker Roth, Eran Halperin, Niko Beerenwinkel
J Comput Biol. 2013 Feb; 20(2): 113–123
