2019 3/18 インストールの流れ修正
2019 3/26 誤字修正
2019 11/8 タイトル修正
ヒトゲノムは二倍体であり、すなわち、その常染色体の各々は2コピーである。これらの親のコピーは、異なる一塩基多型（SNPs）の影響を受ける。変異がどちらの染色体由来かアサインすることは進化遺伝学の助けになり、例えばpopulation研究（論文より The 1000 Genomes Project Consortium、2010; Francioli et al、2014）における選択的な圧力と亜集団の特定に関係し、おそらく病気の原因となるSNPsを相互に結びつける（Hartl and Clark、2007）。アサインするプロセスはphasing（以後フェージング）と呼ばれ、結果のSNPsグループはハプロタイプと呼ばれる。国際的に協調した取り組みにより、対応する下流分析に役立つ多様な集団のハプロタイプの参照パネルが作成された（International HapMap Consortium、2007、2010）。
バリアントのフェージングには2つの主要なアプローチがある。第1のアプローチは、SNPsのアレルのそれらの接合体のステータスつきリストとなる遺伝子型に依存する。ホモ接合の対立遺伝子は相同染色体両コピー上に現れ、明らかに両方のハプロタイプに適用されるが、ヘテロ接合対立遺伝子はコピーのうちの1つにのみ現れ、したがって2つのグループに分けられなければならない。 よってm個のヘテロ接合のSNPsアレルがあると、2のm乗通りのハプロタイプの可能性が生じる。対応する手法は、通常は統計手法になり、既存の参照パネルを統合する。根底にある仮定は、計算されるハプロタイプは、減数分裂中に組換えから生じる基準ハプロタイプブロックのモザイクであるということである。観察された遺伝子型が与えられると、統計的に最も可能性の高いモザイクが出力される。最も一般的なアプローチは、潜在変数モデリング（Howie et al、2009; Li et al、2010; Scheet and Stephens、2006）に基づいている。他のアプローチでは、マルコフ連鎖モンテカルロ法を用いる（Menelaou and Marchini、2013）。
他のアプローチは、シーケンシングデータを直接読み込む。このようなアプローチは、実質的に同一の染色体コピーからのリードを集め、 haplotype assembly approachesと呼ばれる。原則に従えば、目標は、2つのハプロタイプを計算して、最小限の量の配列エラーを訂正してすべてのリードを割り当てることである。このような処方の中で、最近の注目の大半は最小誤差訂正（MEC）問題が生じている。セクション2で正式に定義するMECの問題は、リードを競合のない2つのハプロタイプに配置するために、配列決定されたヌクレオチドに対して行われる補正の最小数を見つけることにある。 MECの主な利点は、phredスケールのエラー確率に対処するために、重み付きバージョン（wMEC）に容易に適合させることができることである。このようなphredベースのエラースキームは、特に将来の技術によって生成されたロングリードを処理するためには不可欠である。 wMEC問題の最適解は、修正されるべきエラーと比較した最尤シナリオに変換される。
ますます増大するリード長および減少するシーケンス決定コストのテラスケールシーケンスプロジェクト（例えば、Boomsma et al、2013; The 1000 Genomes Project Consortium、2010）は、明らかにリードデータから直接フェーズを合わせることが望ましい。しかし、（i）大部分のNGSリードは、いわゆるvariant desertsをブリッジするにはまだ短すぎるため、統計的アプローチが依然として選択される。しかし、連続的なリードベースのフェージングには、隣接するヘテロ接合のSNPsアレルが全てカバーされないといけない。 （ii）MECの問題はNP困難であり、他のすべての同様の問題の定式化もそうである。
MEC（Chen et al、2013; He et al。、2010）の最も進んだ既存のアルゴリズムソリューションは、皮肉なことにvariant desertsから正確に利益を得ていることが多い。しかし、リードベースのアプローチの主な動機は、可能な限り多くのバリアントをカバーするロングリードを処理することで、すべてのバリアントをブリッジすることである。したがって、現在のハプロタイプアセンブリの認識は、それが克服するのが難しい理論上の限界の根底にあることが多い。
この論文では、カバレッジ、つまりSNPs位置をカバーするフラグメントの数が唯一のパラメータである、wMECに対する固定パラメータの扱いやすい（FPT）アプローチであるWhatsHapを紹介する。よって著者らのアプローチの実行時間は、SNPs数の多項式（実際には線形）である。 wMECのための線形ランタイムソリューションは、将来のシークエンシング技術が数万塩基対（bp）のリードを生成し、それらが高い配列決定エラー率を被る可能性が高いことの両方に対処する。慎重に設計されたアルゴリズムを実装することで、標準的なワークステーション上で最大20時間のオーダーで最大カバレッジの全ゲノムデータセットを処理することができる。カバレッジの高いデータセットでは、合理的な選択肢を選択するテクニックを提供する。 WhatsHapはさまざまな面でさまざまな長さのシミュレートされたリードに基づいて比較され、評価された。これらのすべてが真のゲノムに由来する（Levy et al、2007）。そしてハプロタイプアセンブリでのフェージングのエラー率を低くするには10〜15のカバレッジで十分であることを示す。最先端のハプロタイプアセンブリアプローチと比較して、WhatsHapは、特に10〜15カバレッジのデータセット、すなわち重要なカバレッジレートに優れたランタイムを持つことを実証する。 WhatsHapによって計算されたハプロタイプはわずかな誤差しかないことが証明されている。統計的なフェージングアプローチと比較してもエラー率がかなり有利であることは注目に値するかもしれないが、依然として現実的な最先端技術を構成している。
マニュアル
https://whatshap.readthedocs.io/en/latest/guide.html
イルミナ、Nanopore、pacbioいずれのデータセットでも機能する。
インストール
Bitbucket
https://bitbucket.org/whatshap/whatshap
condaによるインストールと、pipによるインストールがサポートされている。
#Anacondaを使っているなら
conda install -c bioconda -c conda-forge whatshap
>whatshap -h
$ whatshap -h
usage: whatshap [-h] [--version] [--debug]
{phase,stats,compare,hapcut2vcf,unphase,haplotag,genotype} ...
positional arguments:
{phase,stats,compare,hapcut2vcf,unphase,haplotag,genotype}
phase Phase variants in a VCF with the WhatsHap algorithm
stats Print phasing statistics of a single VCF file
compare Compare two or more phasings
hapcut2vcf Convert hapCUT output format to VCF
unphase Remove phasing information from a VCF file
haplotag Tag reads by haplotype
genotype Genotype variants
optional arguments:
-h, --help show this help message and exit
--version show program's version number and exit
--debug Print debug messages
> whatshap phase
$ whatshap phase
usage: whatshap phase [-h] [-o OUTPUT] [--reference FASTA] [--tag {PS,HP}]
[--output-read-list FILE]
[--algorithm {whatshap,hapchat}] [--merge-reads]
[--max-coverage MAXCOV] [--mapping-quality QUAL]
[--indels] [--ignore-read-groups] [--sample SAMPLE]
[--chromosome CHROMOSOME]
[--error-rate READ_MERGING_ERROR_RATE]
[--maximum-error-rate READ_MERGING_MAX_ERROR_RATE]
[--threshold READ_MERGING_POSITIVE_THRESHOLD]
[--negative-threshold READ_MERGING_NEGATIVE_THRESHOLD]
[--full-genotyping] [--distrust-genotypes]
[--include-homozygous] [--default-gq DEFAULT_GQ]
[--gl-regularizer GL_REGULARIZER]
[--changed-genotype-list FILE] [--ped PED/FAM]
[--recombination-list FILE] [--recombrate RECOMBRATE]
[--genmap FILE] [--no-genetic-haplotyping]
[--use-ped-samples]
VCF [PHASEINPUT [PHASEINPUT ...]]
Phase variants in a VCF with the WhatsHap algorithm
Read a VCF and one or more files with phase information (BAM/CRAM or VCF phased
blocks) and phase the variants. The phased VCF is written to standard output.
positional arguments:
VCF VCF file with variants to be phased (can be gzip-
compressed)
PHASEINPUT BAM, CRAM or VCF file(s) with phase information,
either through sequencing reads (BAM/CRAM) or through
phased blocks (VCF)
optional arguments:
-h, --help show this help message and exit
-o OUTPUT, --output OUTPUT
Output VCF file. Add .gz to the file name to get
compressed output. If omitted, use standard output.
--reference FASTA, -r FASTA
Reference file. Provide this to detect alleles through
re-alignment. If no index (.fai) exists, it will be
created
--tag {PS,HP} Store phasing information with PS tag (standardized)
or HP tag (used by GATK ReadBackedPhasing) (default:
PS)
--output-read-list FILE
Write reads that have been used for phasing to FILE.
--algorithm {whatshap,hapchat}
Choose an algorithm from whatshap or hapchat (default:
whatshap)
Input pre-processing, selection and filtering:
--merge-reads Merge reads which are likely to come from the same
haplotype (default: do not merge reads)
--max-coverage MAXCOV, -H MAXCOV
Coverage reduction parameter in the core phasing
algorithm. Phasing quality may improve slightly if
this is increased, but runtime increases
exponentially! Do not increase beyond 20. (default:
15)
--mapping-quality QUAL, --mapq QUAL
Minimum mapping quality (default: 20)
--indels Also phase indels (default: do not phase indels)
--ignore-read-groups Ignore read groups in BAM/CRAM header and assume all
reads come from the same sample.
--sample SAMPLE Name of a sample to phase. If not given, all samples
in the input VCF are phased. Can be used multiple
times.
--chromosome CHROMOSOME
Name of chromosome to phase. If not given, all
chromosomes in the input VCF are phased. Can be used
multiple times.
Read merging:
The options in this section are only active when --merge-reads is used
--error-rate READ_MERGING_ERROR_RATE
The probability that a nucleotide is wrong in read
merging model (default: 0.15).
--maximum-error-rate READ_MERGING_MAX_ERROR_RATE
The maximum error rate of any edge of the read merging
graph before discarding it (default: 0.25).
--threshold READ_MERGING_POSITIVE_THRESHOLD
The threshold of the ratio between the probabilities
that a pair of reads come from the same haplotype and
different haplotypes in the read merging model
(default: 1000000).
--negative-threshold READ_MERGING_NEGATIVE_THRESHOLD
The threshold of the ratio between the probabilities
that a pair of reads come from different haplotypes
and the same haplotype in the read merging model
(default: 1000).
Genotyping:
The options in this section require that either --distrust-genotypes or --full-genotyping is used
--full-genotyping Completely re-genotype all variants based on read
data, ignores all genotype data that might be present
in the VCF (EXPERIMENTAL FEATURE).
--distrust-genotypes Allow switching variants from hetero- to homozygous in
an optimal solution (see documentation).
--include-homozygous Also work on homozygous variants, which might be
turned to heterozygous
--default-gq DEFAULT_GQ
Default genotype quality used as cost of changing a
genotype when no genotype likelihoods are available
(default 30)
--gl-regularizer GL_REGULARIZER
Constant (float) to be used to regularize genotype
likelihoods read from input VCF (default None).
--changed-genotype-list FILE
Write list of changed genotypes to FILE.
Pedigree phasing:
--ped PED/FAM Use pedigree information in PED file to improve
phasing (switches to PedMEC algorithm). Columns 2, 3,
4 must refer to child, mother, and father sample names
as used in the VCF and BAM/CRAM. Other columns are
ignored.
--recombination-list FILE
Write putative recombination events to FILE.
--recombrate RECOMBRATE
Recombination rate in cM/Mb (used with --ped). If
given, a constant recombination rate is assumed
(default: 1.26cM/Mb).
--genmap FILE File with genetic map (used with --ped) to be used
instead of constant recombination rate, i.e. overrides
option --recombrate.
--no-genetic-haplotyping
Do not merge blocks that are not connected by reads
(i.e. solely based on genotype status). Default: when
in --ped mode, merge all blocks that contain at least
one homozygous genotype in at least one individual
into one block.
--use-ped-samples Only work on samples mentioned in the provided PED
file.
whatshap phase: error: the following arguments are required: VCF, PHASEINPUT
ラン
whatshapのphaseコマンドは以下のような処理を行う。
これはSample1とSample2のマージしたvcfファイルだが、samle1だけ、GT（genotype）フィールドが５つのアレルすべてでヘテロ（0/1）になっている。
公式より転載
phaseコマンドはGTフィールドが0/1のものをフェーズがそろうよう並べ替えて、/ のかわりにパイプ（ | ）で区切る。Sample1の１つのハプロタイプはA, G, T, T, A、もう片方のハプロタイプはT, C, G, A, Gと推定した（0/1ではヘテロとしかわからない。|で区切ることで、REFかALTか明確にしている）。
公式より転載
phasingを行うには、バリアントコーラーが出力したVCFと、マッピングしたbamファイルが必要になる。 phaseコマンドを使う。
whatshap phase -o phased.vcf input.vcf mapped.bam
1000genomesのデータを３つ解析してみる。サンプル名は通し番号とする。
シェルスクリプト variant_call.sh
#!/bin/sh
#複数のbamを順に処理する。
#依存ツール
#BBtools - interleave fastqの作成
#sickle - quality trimming
#SAMTools - sam=>bam、ソート、index
#BWA MEM - mapping
#vcftools - VCFのマージ
#samblaster
#パラメータ
thread=20 #スレッド数
#保存場所
#全サンプルのfastq: ~/fastq/ #ペアエンド名はstrain通し番号_pair1/2.fq.gz (e.g., strain1_pair1.fq.gz & strain1_pair2.fq.gz)
#interleave_fastq: ~/interleave_fastq/
#reference.fasta: ~/ref/
#bam: ~/bam/
#VCF: ~/vcf/
#indexを作る。
ref=~/ref/Homo_sapiens.GRCh38.dna.chromosome.fa
bwa index $ref #isはたしか2GBまで
samtools faidx $ref
#3サンプルのループ処理
sample=(1 2 3)
for name in ${sample[@]}; do
echo sample$name
r1=strain${name}_pair1.fq.gz
r2=strain${name}_pair2.fq.gz
#interleaveのfastq.gz作成#interleaveのfastq.gz作成
reformat.sh in1=$r1 in2=$r2 out=~/interleave_fastq/${name}_interleave.fq.gz
fastq=~/interleave_fastq/${name}_interleave.fq.gzfastq=~/interleave_fastq/${name}_interleave.fq.gz
sleep 1s
#quality trimming
sickle pe -c $fastq -g -t sanger -l 20 -q 20 -m ~/interleave_fastq/${name}_interleave_qc.fq.gz -s single.fq.gz && rm single.fq.gz
qcfastq=~/interleave_fastq/${name}_interleave_qc.fq.gz
#mapping samblasrにパイプしてduplicationは捨てる
bwa mem -R "@RG\tID:1a\tLB:library1\tSM:Sample${name}\tPL:ILLUMINA" -t $thread $ref $qcfastq | samblaster -r -e -d samp.disc.sam -s samp.split.sam |samtools view -@ $thread -bS - > ${name}.bam
samtools sort -@ $thread ${name}.bam > ~/bam/${name}_sorted.bam && rm ${name}.bam #ソート
samtools index ~/bam/${name}_sorted.bam #bamのindex
freebayes -f $ref -C 5 -p 2 ~/bam/${name}_sorted.bam > ~/vcf/${name}_output.vcf && #バリアントのコール
bgzip -c ~/vcf/${name}_output.vcf > ~/vcf/${name}_output.vcf.gz && rm ~/vcf/${name}_output.vcf #vcfの圧縮
tabix -p vcf ~/vcf/${name}_output.vcf.gz #vcfのindex
done
GATKは遅すぎるのでfreebayesを使った。1_sorted.bam、2_sorted.bam、3_sorted.bam、とそのvcfファイルである1_output.vcf.gz、2_output.vcf.gz、3_output.vcf.gzが得られた。
VCFのマージ。
vcf-merge 1_output.vcf.gz 2_output.vcf.gz 3_output.vcf.gz | bgzip -c > merge.vcf.gz
merge.vcf.gzが出力される。
Sample1のphasingを行う。 WhatsHapのphaseコマンドを使う。ロングリードならリファンレスを与えて精度をあげる必要がある（see manual）。Sample1のchromosome1のみに限定する。
whatshap phase --indels --sample SampleA --chromosome 1 -o phased.vcf merge.vcf.gz 1_sorted.bam
- --max-coverage Reduce coverage (default: 15).
- --mapping-quality Minimum mapping quality (default: 20)
- --indels Also phase indels (default: do not phase indels)
- --ignore-read-groups Ignore read groups in BAM header and assume all reads come from the same sample.
- --sample Name of a sample to phase. If not given, all samples in the input VCF are phased. Can be used multiple
- --chromosome Name of chromosome to phase. If not given, all chromosomes in the input VCF are phased. Can be used multiple times.
phased.vcfが出力される。
gz圧縮してindexをつける。
bgzip phased.vcf
tabix phased.vcf.gz
フェーズ分析結果をもとにハプロタイプの配列を出力。
bcftools consensus -H 1 -s Sample1 -f Homo_sapiens.GRCh38.dna.chromosome.fa phased.vcf.gz > haplotype1.fa
bcftools consensus -H 2 -s Sample1 -f Homo_sapiens.GRCh38.dna.chromosome.fa phased.vcf.gz > haplotype2.fa
フェーズ分析結果をもとにGTFを出力。phasingされた数などのstatisticsも標準出力される。
whatshap stats --sample SampleA --chromosome 1 --indels --gtf=phased.gtf phased.vcf.gz
- --gtf Write phased blocks to GTF file.
- --sample Name of the sample to process. If not given, use first sample found in VCF.
- --only-snvs Only process SNVs and ignore all other variants.
- --chromosome Name of chromosome to process. If not given, all chromosomes in the input VCF are considered.
IGVに読み込む。
igv -g Homo_sapiens.GRCh38.dna.chromosome.fa phased.vcf.gz,phased.gtf,1_sorted.bam
他にも２つ以上のphasing解析結果を比較するcompareなどのコマンドがある。
引用
WhatsHap: Weighted Haplotype Assembly for Future-Generation Sequencing Reads.
Patterson M, Marschall T, Pisanti N, van Iersel L, Stougie L, Klau GW, Schönhuth A
J Comput Biol. 2015 Jun;22(6):498-509.
