#!/bin/bash

set -euo pipefail

ln -sf ../11.reference_v2/raw .
cd raw
ln -sf Hippophae_tibetana.genomic.fasta HT.fa
ln -sf YN.genome.fa YN.fa
seqkit grep -f ZG.target.txt ZG.genome.fa >  ZG.fa
ln -sf ZJ.genome.fa ZJ.fa
seqkit grep -f MG.target.txt MG.genome.fasta >  MG.fa
ln -sf EM.genome.fasta EM.fa
cd ..

bash "$SCRIPTS"/utils/asm_qc.sh raw/ZJ.fa \
  raw/EM.fa \
  raw/EA.fa \
  raw/HT.fa \
  raw/ZY.fa \
  raw/ZG.fa \
  raw/MG.fa \
  raw/YN.fa \
  raw/HG.fa \
  raw/HS.fa \
  raw/HN.fa > qc.tsvWDw

##
## minimap2 比对：YN 为参考（ZG 含单倍型出局），asm5 档位（属内分歧 0.001-0.007，EA/EM 约 0.03 也覆盖）
## ZJ 已剔除（跨科，分歧 0.19）；YN 不做自比对
##
REF=raw/YN.fa
mkdir -p paf
pueue add -- minimap2 -cx asm5 --cs=long -t 8 "$REF" raw/HG.fa '>' paf/HG.paf
pueue add -- minimap2 -cx asm5 --cs=long -t 8 "$REF" raw/HN.fa '>' paf/HN.paf
pueue add -- minimap2 -cx asm5 --cs=long -t 8 "$REF" raw/HS.fa '>' paf/HS.paf
pueue add -- minimap2 -cx asm5 --cs=long -t 8 "$REF" raw/HT.fa '>' paf/HT.paf
pueue add -- minimap2 -cx asm5 --cs=long -t 8 "$REF" raw/MG.fa '>' paf/MG.paf
pueue add -- minimap2 -cx asm5 --cs=long -t 8 "$REF" raw/ZG.fa '>' paf/ZG.paf
pueue add -- minimap2 -cx asm5 --cs=long -t 8 "$REF" raw/ZY.fa '>' paf/ZY.paf
pueue add -- minimap2 -cx asm5 --cs=long -t 8 "$REF" raw/EA.fa '>' paf/EA.paf
pueue add -- minimap2 -cx asm5 --cs=long -t 8 "$REF" raw/EM.fa '>' paf/EM.paf

## 覆盖速览（pilot 判档用；asm5 若 EA/EM 覆盖塌，对这两个换 asm10 重跑）
for f in paf/*.paf; do
    awk -v OFS='\t' '{alen += $11; qlen = $2} END {printf "%s\tqlen=%d\taligned=%d\tcov=%.2f%%\n", FILENAME, qlen, alen, alen / qlen * 100}' "$f"
done

## YN 自比对：Dsuite 需要 YN 作样本（自比对产生全纯合参考基因型 VCF）
pueue add -- minimap2 -cx asm5 --cs=long -t 8 "$REF" raw/YN.fa '>' paf/YN.paf

##
## PAF 过滤：只留 primary 链（tp:A:P）并按参考坐标排序（paftools call 要求排序输入）
## primary 链是一对一共线块的近似；如需严格双向最佳命中，在此之上再加过滤
##
mkdir -p paf_primary
grep 'tp:A:P' paf/HG.paf | sort -k6,6 -k8,8n > paf_primary/HG.paf
grep 'tp:A:P' paf/HN.paf | sort -k6,6 -k8,8n > paf_primary/HN.paf
grep 'tp:A:P' paf/HS.paf | sort -k6,6 -k8,8n > paf_primary/HS.paf
grep 'tp:A:P' paf/HT.paf | sort -k6,6 -k8,8n > paf_primary/HT.paf
grep 'tp:A:P' paf/MG.paf | sort -k6,6 -k8,8n > paf_primary/MG.paf
grep 'tp:A:P' paf/ZG.paf | sort -k6,6 -k8,8n > paf_primary/ZG.paf
grep 'tp:A:P' paf/ZY.paf | sort -k6,6 -k8,8n > paf_primary/ZY.paf
grep 'tp:A:P' paf/YN.paf | sort -k6,6 -k8,8n > paf_primary/YN.paf
grep 'tp:A:P' paf/EA.paf | sort -k6,6 -k8,8n > paf_primary/EA.paf
grep 'tp:A:P' paf/EM.paf | sort -k6,6 -k8,8n > paf_primary/EM.paf

##
## paftools call：每样本一个 VCF（-s 样本名；-l/-L 最小块长阈值，可调）
##
mkdir -p vcf
paftools.js call -s HG -l 5000 -L 10000 paf_primary/HG.paf > vcf/HG.vcf
paftools.js call -s HN -l 5000 -L 10000 paf_primary/HN.paf > vcf/HN.vcf
paftools.js call -s HS -l 5000 -L 10000 paf_primary/HS.paf > vcf/HS.vcf
paftools.js call -s HT -l 5000 -L 10000 paf_primary/HT.paf > vcf/HT.vcf
paftools.js call -s MG -l 5000 -L 10000 paf_primary/MG.paf > vcf/MG.vcf
paftools.js call -s ZG -l 5000 -L 10000 paf_primary/ZG.paf > vcf/ZG.vcf
paftools.js call -s ZY -l 5000 -L 10000 paf_primary/ZY.paf > vcf/ZY.vcf
paftools.js call -s YN -l 5000 -L 10000 paf_primary/YN.paf > vcf/YN.vcf
paftools.js call -s EA -l 5000 -L 10000 paf_primary/EA.paf > vcf/EA.vcf
paftools.js call -s EM -l 5000 -L 10000 paf_primary/EM.paf > vcf/EM.vcf

bgzip -f vcf/HG.vcf && tabix -f -p vcf vcf/HG.vcf.gz
bgzip -f vcf/HN.vcf && tabix -f -p vcf vcf/HN.vcf.gz
bgzip -f vcf/HS.vcf && tabix -f -p vcf vcf/HS.vcf.gz
bgzip -f vcf/HT.vcf && tabix -f -p vcf vcf/HT.vcf.gz
bgzip -f vcf/MG.vcf && tabix -f -p vcf vcf/MG.vcf.gz
bgzip -f vcf/ZG.vcf && tabix -f -p vcf vcf/ZG.vcf.gz
bgzip -f vcf/ZY.vcf && tabix -f -p vcf vcf/ZY.vcf.gz
bgzip -f vcf/YN.vcf && tabix -f -p vcf vcf/YN.vcf.gz
bgzip -f vcf/EA.vcf && tabix -f -p vcf vcf/EA.vcf.gz
bgzip -f vcf/EM.vcf && tabix -f -p vcf vcf/EM.vcf.gz

##
## 合并为 10 类群多样本 VCF
##
pueue add -- bcftools merge -O z -o vcf/dsuite.10taxa.vcf.gz vcf/HG.vcf.gz vcf/HN.vcf.gz vcf/HS.vcf.gz vcf/HT.vcf.gz vcf/MG.vcf.gz vcf/ZG.vcf.gz vcf/ZY.vcf.gz vcf/YN.vcf.gz vcf/EA.vcf.gz vcf/EM.vcf.gz
tabix -f -p vcf vcf/dsuite.10taxa.vcf.gz

##
## 物种树：主数据集第二套（exclude_all_fs）ASTRAL 树剔除 ZJ
## TODO: 剪树工具任选（nw_prune / gotree / R ape），产物放 tree/dsuite.nwk
##
mkdir -p tree
# nw_prune ../14.taxa_phylogeny/exclude_all_fs/species_tree/coal.tre ZJ > tree/dsuite.nwk

##
## SETS 文件：双外类群交叉验证，EA / EM 各标一次 Outgroup
##
mkdir -p sets dsuite
cat > sets/EA_out.txt <<'EOF'
HG	HG
HN	HN
HS	HS
HT	HT
MG	MG
ZG	ZG
ZY	ZY
YN	YN
EA	Outgroup
EM	EM
EOF
cat > sets/EM_out.txt <<'EOF'
HG	HG
HN	HN
HS	HS
HT	HT
MG	MG
ZG	ZG
ZY	ZY
YN	YN
EA	EA
EM	Outgroup
EOF

##
## Dtrios + Fbranch（EA / EM 各跑一遍）
##
pueue add -- Dsuite Dtrios -t tree/dsuite.nwk -o dsuite/EA_out vcf/dsuite.10taxa.vcf.gz sets/EA_out.txt
pueue add -- Dsuite Dtrios -t tree/dsuite.nwk -o dsuite/EM_out vcf/dsuite.10taxa.vcf.gz sets/EM_out.txt
Dsuite Fbranch tree/dsuite.nwk dsuite/EA_out_tree.txt > dsuite/EA_out_fbranch.txt
Dsuite Fbranch tree/dsuite.nwk dsuite/EM_out_tree.txt > dsuite/EM_out_fbranch.txt
# dtools.py dsuite/EA_out_fbranch.txt tree/dsuite.nwk
# dtools.py dsuite/EM_out_fbranch.txt tree/dsuite.nwk