长读测序单倍型组装 phasing¶

# ===== 长读测序 Phasing 流程 =====

# ===== 方法1：WhatsHap — 读基础 phasing =====

# 1. 先用 DeepVariant 等工具 call 变异（SNP + Indel）
# 假设已有 variants.vcf.gz 和 aligned.sorted.bam

# 2. WhatsHap phase：根据长读连接杂合位点
whatshap phase \
    --reference hg38.fa \
    --output phased.vcf.gz \
    variants.vcf.gz \
    aligned.sorted.bam

# 3. 查看 phasing 统计
whatshap stats phased.vcf.gz \
    --gtf gencode.annotation.gtf \
    > phasing_stats.txt

# 关键指标：
# - Phase block N50: 最大连续 phasing 区域的 N50
# - Phased variants: 被 phased 的变异比例
# - Switch error rate: phasing 错误率

# 4. 用 phased VCF 标记 BAM（haplotag）
# 给每条 read 打上 HP:1 或 HP:2 标签
whatshap haplotag \
    --reference hg38.fa \
    --output haplotagged.bam \
    phased.vcf.gz \
    aligned.sorted.bam
samtools index haplotagged.bam

# 5. 按单倍型分割 BAM
# 分出 hap1 和 hap2 的 reads
samtools view -h -d HP:1 haplotagged.bam | \
    samtools sort -@ 4 -o hap1.bam
samtools view -h -d HP:2 haplotagged.bam | \
    samtools sort -@ 4 -o hap2.bam
samtools index hap1.bam
samtools index hap2.bam

# ===== 方法2：hifiasm 直接单倍型组装 =====

# 使用 Hi-C 数据辅助 phasing
hifiasm -o asm_phased -t 32 \
    --h1 hic_R1.fastq.gz \
    --h2 hic_R2.fastq.gz \
    hifi_reads.fastq.gz

# 输出两套单倍型组装
awk '/^S/{print ">"$2; print $3}' \
    asm_phased.hic.hap1.p_ctg.gfa > hap1_assembly.fa
awk '/^S/{print ">"$2; print $3}' \
    asm_phased.hic.hap2.p_ctg.gfa > hap2_assembly.fa

# ===== 方法3：Trio phasing（需要父母数据） =====

# 1. 用 yak 建父母 k-mer 库
yak count -t 16 -b 37 -o pat.yak father_illumina.fastq.gz
yak count -t 16 -b 37 -o mat.yak mother_illumina.fastq.gz

# 2. hifiasm trio 模式
hifiasm -o asm_trio -t 32 \
    -1 pat.yak -2 mat.yak \
    child_hifi.fastq.gz

# 输出：完全按父母来源分开的两套组装
# asm_trio.dip.hap1.p_ctg.gfa → 父源单倍型
# asm_trio.dip.hap2.p_ctg.gfa → 母源单倍型

# ===== Python: Phasing 质量评估 =====
import pysam
from collections import Counter

# 统计 haplotagged BAM 的 phasing 情况
bam = pysam.AlignmentFile("haplotagged.bam", "rb")

hp_counts = Counter()
unphased = 0
total = 0

for read in bam:
    if read.is_unmapped or read.is_secondary:
        continue
    total += 1
    if read.has_tag("HP"):
        hp = read.get_tag("HP")  # 1 或 2
        hp_counts[hp] += 1
    else:
        unphased += 1

bam.close()

print(f"总 reads 数: {total}")
print(f"Hap1 reads: {hp_counts[1]} ({hp_counts[1]/total*100:.1f}%)")
print(f"Hap2 reads: {hp_counts[2]} ({hp_counts[2]/total*100:.1f}%)")
print(f"未 phased: {unphased} ({unphased/total*100:.1f}%)")
# 理想情况：Hap1 和 Hap2 大致相等，未 phased 比例低

# 解析 phasing 统计文件
print("\n=== Phase Block 统计 ===")
with open("phasing_stats.txt") as f:
    for line in f:
        if "N50" in line or "phased" in line.lower():
            print(line.strip())