189_全基因组倍性分析¶

# 计算K-mer频率
jellyfish count -C -m 21 -s 10G -t 32 reads_R1.fq.gz reads_R2.fq.gz -o counts.jf
jellyfish histo -t 8 counts.jf > kmer_histo.txt

# 用GenomeScope2自动尝试不同倍性
for p in 1 2 3 4 6; do
    genomescope2 \
        -i kmer_histo.txt \
        -o genomescope_p${p} \
        -k 21 \
        -p ${p} \
        --testing
done

# 比较不同倍性模型的拟合优度
for p in 1 2 3 4 6; do
    echo "Ploidy ${p}:"
    grep "GenomeScope" genomescope_p${p}/summary.txt
done

# Smudgeplot验证
# 提取杂合K-mer对
jellyfish dump -c -L 2 counts.jf > kmer_dump.txt
smudgeplot.py hetkmers -o pairs kmer_dump.txt
smudgeplot.py plot pairs_coverages.tsv -o smudgeplot_result

# Smudgeplot会显示最可能的倍性和预测置信度

# nQuire基于BAF分析
nQuire create -b sample.bam -o sample
nQuire denoise sample.bin -o sample_dn
nQuire lrdmodel sample_dn.bin

# ploidyNGS
ploidyNGS --bam sample.bam --out ploidy_result

# ConPADE（基于等位基因计数的倍性估计）
python conpade.py -b sample.bam -r reference.fa -o conpade_result

# 流式数据分析（R）
library(flowCore)

# 读取FCS文件
fcs <- read.FCS("ploidy_analysis.fcs")

# 提取PI/DAPI荧光通道
dna_signal <- exprs(fcs)[, "FL2-A"]  # 根据具体通道调整

# 绘制直方图
hist(dna_signal, breaks = 200, main = "DNA Content Distribution",
     xlab = "Fluorescence Intensity", col = "steelblue", border = NA)

# 与对照比较
# 如果样品的G1峰是对照（2n）的1.5倍 → 三倍体
# 如果样品的G1峰是对照（2n）的2倍 → 四倍体

import matplotlib.pyplot as plt
import numpy as np

# 不同倍性的K-mer频谱模拟
fig, axes = plt.subplots(1, 4, figsize=(20, 4))
x = np.linspace(0, 200, 1000)
cov = 50

titles = ['Haploid (1n)', 'Diploid (2n)', 'Triploid (3n)', 'Tetraploid (4n)']
for idx, (ploidy, ax, title) in enumerate(zip([1, 2, 3, 4], axes, titles)):
    y = np.zeros_like(x)
    for k in range(1, ploidy + 1):
        peak = cov * k / ploidy
        y += np.exp(-(x - peak)**2 / (2 * (peak * 0.15)**2)) * (ploidy - k + 1)
    ax.plot(x, y, 'b-', linewidth=2)
    ax.fill_between(x, y, alpha=0.3)
    ax.set_title(title)
    ax.set_xlabel('K-mer Coverage')
    ax.set_ylabel('Frequency')
    for k in range(1, ploidy + 1):
        peak = cov * k / ploidy
        ax.axvline(x=peak, color='red', linestyle='--', alpha=0.5)

plt.tight_layout()
plt.savefig('ploidy_kmer_patterns.png', dpi=300)

# 完整倍性检测流程
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 2  # 测试二倍体
genomescope2 -i histo.txt -o gs_out3 -k 21 -p 3  # 测试三倍体
smudgeplot.py hetkmers -o pairs kmer_dump.txt
smudgeplot.py plot pairs_coverages.tsv -o smudge

# 基于BAM的快速检测
nQuire create -b aligned.bam -o sample; nQuire lrdmodel sample.bin