178_全基因组三倍体分析¶

# 1. 计算k-mer频率（使用jellyfish）
jellyfish count -C -m 21 -s 5G -t 16 -o kmer_counts.jf reads_R1.fq.gz reads_R2.fq.gz
jellyfish histo -t 16 kmer_counts.jf > kmer_histo.txt

# 2. 使用GenomeScope2估计倍性和基因组大小
# 安装：https://github.com/tbenavi1/genomescope2.0
genomescope2 -i kmer_histo.txt -o genomescope_output -k 21 -p 3
# -p 3 指定三倍体模型

# 也可以自动尝试不同倍性
for ploidy in 2 3 4 6; do
    genomescope2 -i kmer_histo.txt -o genomescope_p${ploidy} -k 21 -p ${ploidy}
done

# 3. 使用Smudgeplot进一步验证倍性
# Smudgeplot通过k-mer对的频率比来估计倍性
# 提取k-mer pairs
jellyfish count -C -m 21 -s 5G -t 16 -o kmer_21.jf reads_*.fq.gz
jellyfish dump -c -L 2 kmer_21.jf | sort -k2,2nr > kmer_dump.txt

# 计算smudgeplot
smudgeplot.py hetkmers -o kmer_pairs kmer_dump.txt
smudgeplot.py plot kmer_pairs_coverages.tsv -o smudgeplot -t "Sample Smudgeplot"

# 安装nQuire
git clone https://github.com/clwgg/nQuire.git
cd nQuire && make

# 1. 将BAM文件转换为nQuire输入
./nQuire create -b sample.bam -o sample

# 2. 去噪
./nQuire denoise sample.bin -o sample_denoised

# 3. 运行倍性估计
./nQuire lrdmodel sample_denoised.bin
# 输出各倍性模型的delta log-likelihood
# diploid, triploid, tetraploid

# 4. 查看BAF分布
./nQuire histo sample_denoised.bin > baf_histo.txt

# 可视化BAF分布
python3 << 'EOF'
import matplotlib.pyplot as plt
import numpy as np

data = np.loadtxt("baf_histo.txt")
plt.figure(figsize=(8, 5))
plt.bar(data[:, 0], data[:, 1], width=0.01, alpha=0.7)
plt.xlabel("B-Allele Frequency")
plt.ylabel("Count")
plt.title("BAF Distribution")
# 三倍体应该在0.33和0.67处有峰
plt.axvline(x=0.33, color='r', linestyle='--', label='1/3')
plt.axvline(x=0.67, color='r', linestyle='--', label='2/3')
plt.legend()
plt.savefig("baf_distribution.png", dpi=300)
EOF

# 使用hifiasm组装三倍体基因组
hifiasm -o output \
    -t 32 \
    --n-hap 3 \
    hifi_reads.fq.gz

# 如果有Hi-C数据辅助分型
hifiasm -o output \
    -t 32 \
    --n-hap 3 \
    --h1 hic_R1.fq.gz \
    --h2 hic_R2.fq.gz \
    hifi_reads.fq.gz

# 提取contigs
awk '/^S/{print ">"$2; print $3}' output.p_ctg.gfa > primary.fa
awk '/^S/{print ">"$2; print $3}' output.hap1.p_ctg.gfa > hap1.fa
awk '/^S/{print ">"$2; print $3}' output.hap2.p_ctg.gfa > hap2.fa

# 使用purge_dups去冗余
minimap2 -xmap-hifi primary.fa hifi_reads.fq.gz | gzip > aligned.paf.gz
pbcstat aligned.paf.gz
calcuts PB.stat > cutoffs
purge_dups -2 -T cutoffs -c PB.base.cov primary.fa > purged.fa

# 评估组装质量
quast purged.fa -o quast_output -t 16
busco -i purged.fa -l embryophyta_odb10 -m genome -c 16

# GATK HaplotypeCaller设置三倍体
gatk HaplotypeCaller \
    -R reference.fa \
    -I sample.bam \
    -O sample.g.vcf.gz \
    -ERC GVCF \
    --sample-ploidy 3

# GenotypeGVCFs
gatk GenotypeGVCFs \
    -R reference.fa \
    -V sample.g.vcf.gz \
    -O genotyped.vcf.gz \
    --sample-ploidy 3

# 过滤变异
gatk VariantFiltration \
    -R reference.fa \
    -V genotyped.vcf.gz \
    -O filtered.vcf.gz \
    --filter-expression "QD < 2.0" --filter-name "LowQD" \
    --filter-expression "FS > 60.0" --filter-name "HighFS"

# 查看三倍体基因型格式
# 三倍体SNP示例：0/0/1 (AAB), 0/1/1 (ABB), 0/1/2 (ABC)
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT[\t%GT]\n' filtered.vcf.gz | head

# 完整倍性鉴定流程
jellyfish count -C -m 21 -s 5G -t 16 reads.fq.gz -o counts.jf
jellyfish histo counts.jf > histo.txt
genomescope2 -i histo.txt -o gs_out -k 21 -p 3

# nQuire快速倍性检测
nQuire create -b sample.bam -o sample
nQuire denoise sample.bin -o sample_dn
nQuire lrdmodel sample_dn.bin

# ploidyNGS基于BAM的倍性估计
ploidyNGS --bam sample.bam --out ploidy_result

# 三倍体GATK流程
gatk HaplotypeCaller -R ref.fa -I sample.bam -O out.g.vcf.gz -ERC GVCF --sample-ploidy 3