143 lines
5.4 KiB
Bash
143 lines
5.4 KiB
Bash
#!/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 |