基因组Survey分析¶

原始reads → 质控(fastp) → k-mer计数(Jellyfish) → 频率分布 → GenomeScope估算
                                                              ↓
                                                  基因组大小、杂合度、重复率

# === 用fastp做质控 ===
fastp \
    -i raw_R1.fastq.gz \         # 原始正向reads
    -I raw_R2.fastq.gz \         # 原始反向reads
    -o clean_R1.fastq.gz \       # 质控后正向reads
    -O clean_R2.fastq.gz \       # 质控后反向reads
    -q 20 \                      # 最低碱基质量
    -l 50 \                      # 最短读段长度
    -w 8 \                       # 工作线程
    -h fastp_report.html          # HTML报告

# === 使用Jellyfish计数k-mer ===

# 第1步：计数（k=21最常用）
jellyfish count \
    -C \                         # 正义链和反义链合并计数
    -m 21 \                      # k-mer长度=21
    -s 5G \                      # 哈希表初始大小（5GB，大基因组用更大值）
    -t 16 \                      # 线程数
    -o kmer_counts.jf \          # 输出文件
    <(zcat clean_R1.fastq.gz) \  # 解压输入（正向）
    <(zcat clean_R2.fastq.gz)    # 解压输入（反向）

# 第2步：生成频率直方图
jellyfish histo \
    -t 16 \                      # 线程数
    -h 1000000 \                 # 最大频率上限
    kmer_counts.jf > kmer_histo.txt   # 输出直方图

# === 使用KMC替代方案（更快，内存更少） ===
# 适合大基因组（>1Gb）

# 创建临时目录
mkdir kmc_tmp

# 运行KMC
kmc \
    -k21 \                       # k-mer长度
    -ci2 \                       # 最低计数阈值（过滤测序错误）
    -cs1000000 \                 # 最大计数值
    -t16 \                       # 线程数
    -m64 \                       # 最大内存(GB)
    @input_list.txt \            # 输入文件列表（每行一个fastq路径）
    kmc_output \                 # 输出前缀
    kmc_tmp                       # 临时目录

# 生成直方图
kmc_tools transform kmc_output histogram kmer_histo.txt

# === 方法1: GenomeScope 2.0 在线版 ===
# 访问: http://qb.cshl.edu/genomescope/genomescope2.0/
# 上传kmer_histo.txt文件
# 设置k-mer length=21, Ploidy=2(二倍体)
# 点击Submit

# === 方法2: GenomeScope 2.0 命令行 ===
# 安装
conda install -c bioconda genomescope2

# 运行
genomescope2 \
    -i kmer_histo.txt \          # k-mer直方图文件
    -k 21 \                      # k-mer长度
    -p 2 \                       # 倍性（2=二倍体，4=四倍体）
    -o genomescope_output         # 输出目录

# 输出文件:
# genomescope_output/
# ├── summary.txt         # ★文字摘要（基因组大小、杂合度等）
# ├── model.txt           # 模型参数
# ├── linear_plot.png     # ★线性坐标的k-mer分布图
# ├── log_plot.png        # ★对数坐标的k-mer分布图
# └── transformed_linear_plot.png  # 变换后的图

# === GenomeScope summary.txt 示例 ===
# Genome Haploid Length   485,234,567 bp     # 单倍体基因组大小 ★
# Genome Repeat Length    198,765,432 bp     # 重复序列长度
# Genome Unique Length    286,469,135 bp     # 唯一序列长度
# Heterozygosity           1.23 %            # 杂合度 ★
# Repeat                   40.97 %           # 重复序列比例 ★
# Read Error Rate           0.35 %           # 读段错误率
# Model Fit                 97.5 %           # 模型拟合度

频率
 ^
 |    *
 |   * *         ← 杂合峰（在主峰的1/2位置，杂合度高时明显）
 |  *   *
 | *     *****
 |*         * *  ← 主峰（代表单倍体覆盖度）
 |           * *
 |            * *** ← 右侧长尾（重复序列）
 +--+--+--+--+--→ 覆盖度(depth)
    0  25 50 75 100

# 主峰位置 = 平均测序深度
# 杂合峰 = 主峰的一半位置出现的次峰
# 右侧长尾 = 重复序列（出现频率是主峰的2倍、3倍...）
# 最左侧的极高频点 = 测序错误产生的k-mer（应过滤）

# === 基本公式 ===
# 基因组大小 = k-mer总数 / 主峰覆盖度

# 从直方图计算k-mer总数
awk '{total += $1 * $2} END{print "Total k-mers:", total}' kmer_histo.txt

# 假设主峰在覆盖度50处
# 基因组大小 = k-mer总数 / 50

# === 更精确的计算（排除错误k-mer） ===
# 通常覆盖度<5的k-mer被认为是测序错误
awk '$1 >= 5 {total += $1 * $2} END{print "Filtered k-mers:", total}' kmer_histo.txt
# 基因组大小 = 过滤后k-mer总数 / 主峰覆盖度

# === Smudgeplot：可视化基因组倍性和杂合度 ===
# 特别适合判断多倍体和杂合体

# 安装
conda install -c bioconda smudgeplot

# 从KMC结果提取k-mer对
smudgeplot.py hetkmers \
    -o kmer_pairs \              # 输出前缀
    kmc_output                    # KMC数据库

# 绘制smudgeplot
smudgeplot.py plot \
    kmer_pairs_coverages.tsv \   # 覆盖度文件
    -o smudge_plot \             # 输出前缀
    -k 21                        # k-mer长度

# 结果解读：
# AB类型的点集中在 → 二倍体(AB)
# 多个点簇 → 可能是多倍体
# 点分散 → 高杂合度

# === 基因组Survey标准流程 ===
# 1. 质控
fastp -i R1.fq.gz -I R2.fq.gz -o clean_R1.fq.gz -O clean_R2.fq.gz

# 2. k-mer计数
jellyfish count -C -m 21 -s 5G -t 16 -o counts.jf <(zcat clean_*.fq.gz)
jellyfish histo -t 16 counts.jf > histo.txt

# 3. GenomeScope估算
genomescope2 -i histo.txt -k 21 -p 2 -o results/

# 关键输出：基因组大小、杂合度、重复率
# 推荐测序深度：Survey用30x，正式组装用50-100x（NGS）或30x（HiFi）