生信软件37 - 基于测序reads的变异进行单倍型分型工具WhatsHap

1. WhatsHap简介

WhatsHap是一种使用DNA测序reads的基因组变异进行定相（分型）的软件，即基于reads的定相或单倍型组装，特别适用于长reads （三代测序数据），但也兼容短reads的定相。

Whatshap特点：

定相结果准确；
适用于Illumina、PacBio、Oxford Nanopor等测序reads数据；
支持对SNV、indel甚至复杂变异（如TCG → AGAA）进行定相；
谱系定相模式使用来自相关个体（例如家系三人组）的reads来改善结果并降低覆盖要求；
较容易安装；
易于使用：输入1个VCF文件和一个BAM文件或多个BAM文件，输出一个定相的VCF；

2. Whatshap安装

支持conda安装和python安装，但python安装不能安装依赖，推荐conda安装。

# 直接conda安装
conda install whatshap -y


# conda 创建新环境并安装whashap
conda create -n whatshap-env whatshap

# 激活
conda activate whatshap-env

# 查看版本
whatshap --version

# 获取子命令帮助信息
whatshap SUBCOMMAND --help

Whatshap子命令：

![[子命令]]

3. 基于reads的定相和用WhatsHap分析定相变异

文献： https://doi.org/10.1007/978-1-0716-2819-5_8

在线文档：https://whatshap.readthedocs.io/

github地址： https://github.com/whatshap/whatshap.git

![[文献]]

4. 定相基本用法

建议使用高质量的短reads进行变异识别，然后使用长reads进行定相。WhatsHap可以对SNVs（单核苷酸变异）、插入、缺失、MNP（多个相邻SNVs）和“复杂”变异进行定相。软件将输入VCF中标记为杂合（基因型0/1）且具有适当覆盖度的变异作为核心定相算法的输入信息。

如果输入VCF是多样本VCF，WhatsHap将单独对所有样本进行单倍型分型。如果要对家系相关个体的样本进行定相，可以使用谱系定相模式来改进结果。

建议使用--reference提供FASTA参考基因组序列文件，特别是对于易出错的reads（PacBio，Nanopore），以便软件能实现重新比对变异检测算法。如果未提供参考序列文件，则仅可对SNV、插入和缺失进行定相。

# -o: 输出定相vcf文件
# --reference: 参考基因组路径
# input.vcf: 输入相同个体的vcf文件
# input.bam: 输入相同个体的bam文件

whatshap phase \
-o phased.vcf \
--reference=reference.fasta \
input.vcf \
input.bam

5. 谱系定相

Whatshap只使用ped文件的第2列（个体ID）、第3列（父亲ID）、第4列（母亲ID）。
谱系模式需要使用获取重组事件的成本，默认情况下，WhatsHap将假设要定相的染色体上的重组率恒定，可以通过提供选项--recombrate来改变重组率（以cM/Mb为单位）。默认值1.26 cM/Mb适用于人类基因组。

# --ped: PINK格式的ped文件, 至少有六列，#开头内容被忽略
# 示例如下：

# Fields: family, individual_id, paternal_id, maternal_id, sex, phenotype
FAMILY01 child father mother 0 1

whatshap phase \
--ped pedigree.ped \
--reference=reference.fasta \
-o phased.vcf \
input.vcf \
input.bam

6. 统计相位信息

# 打印单个VCF文件的定相统计信息
whatshap stats input.vcf

7. 可视化分型结果

借助IGV 基因组浏览器可视化定相结果。WhatsHap可以从描述单倍型块的定相VCF文件创建GTF文件，只需要使用phased.vcf中的定相结果作为输入。使用打开IGV中的phased.vcf和phased.gtf以检查单倍型块。

whatshap stats \
--gtf=phased.gtf \
phased.vcf

如果配对或双端PE读段用于定相，则单体型块可以是交错的或嵌套的。

![[可视化]]

8. 通过单体型标记读数以实现可视化

# -o: 输出单倍型标记BAM文件
# 在BAM文件中根据其属于哪种单倍型用`HP:i:1`或`HP:i:2`标记每个reads

whatshap haplotag \
-o haplotagged.bam \
--reference reference.fasta \
phased.vcf.gz \
alignments.bam

8.1 可视化单倍型块

IGV右键单击BAM轨迹并选择Color Alignments by → tag, 然后键入PS并单击“确定”。

![[haplotagged-PS.webp]]

8.2 可视化单倍型分配

IGV选择Color Alignments by → tag并键入HP。

红色的读数属于一种单倍型，而蓝色的读数属于另一种单倍型。灰色读数是那些无法标记的读数，通常是因为它们不覆盖任何杂合变异。

![[haplotagged-HP.webp]]

9. 对变异进行基因分型

给定包含变异位置的VCF文件，它计算所有三种基因型（0/0，0/1，1/1）的基因型可能性，并将它们与基因型预测一起输出到VCF文件中。

定相时强烈建议提供参考序列，以启用重新对齐模式。

# -o： 输出包含分型结果的vcf文件

whatshap genotype \
--reference ref.fasta \
-o genotyped.vcf \
variants.vcf \
reads.bam

如果没有输入VCF文件可用，WhatsHap可以产生候选SNV位置，这些位置可以用作上述基因分型命令的输入。以下命令来实现：

whatshap find_snv_candidates \
ref.fasta \
input.bam \
-o variants.vcf


# 如果Nanopore读取用于调用SNP，需将选项-nanopore添加到上述命令中
whatshap find_snv_candidates \
ref.fasta \
input.bam \
-o variants.vcf \
-nanopore

10. 二倍体定相

除了二倍体定相，WhatsHap还通过不同的算法支持多倍体定相。命令与phase模式命令类似。

# --ploidy: 倍性，人类为2倍体

whatshap polyphase \
input.vcf \
input.bam \
--ploidy p \
--reference ref.fasta \
-o output.vcf

11. 比对变异VCF文件

whatshap compare主要根据切换误差（switch errors）来评估差异，但它也计算翻转误差和汉明距离（flip errors and Hamming distance）。

# --names用于将名称“truth”分配给第一个输入文件，将“whatshap”分配给第二个输入文件。
# --tsv-pairwise: \t分割的结果文件

whatshap compare \
--names truth,whatshap \
--tsv-pairwise eval.tsv \
truth.chr1.vcf \
phased.chr1.vcf

12. 使用haplotaged BAM文件对VCF文件进行分型

给定一个单倍标记的BAM文件（haplotagged.bam）、一个VCF文件和一个bam index，此命令序列输出一个分型的VCF文件。

# 常见vcf.gz索引
tabix input.vcf.gz

samtools index haplotagged.bam

# -o: 输出分型的vcf文件

whatshap haplotagphase \
-r reference.fasta \
input.vcf.gz \
haplotagged.bam \
-o output.vcf.gz