基因组从头组装策略¶

# === 使用MEGAHIT进行短读段组装 ===
# 适合：宏基因组、小基因组
megahit \
    -1 reads_R1.fastq.gz \       # 双端测序的正向reads
    -2 reads_R2.fastq.gz \       # 双端测序的反向reads
    --min-contig-len 500 \       # 最短contig长度，过滤掉太短的片段
    --k-list 21,41,61,81,99 \    # 多个k-mer值，从小到大逐步组装
    -t 16 \                      # 使用16个CPU线程
    -o megahit_output             # 输出目录

# === 使用SPAdes进行短读段组装 ===
# 适合：细菌基因组、小真核基因组
spades.py \
    -1 reads_R1.fastq.gz \       # 正向reads
    -2 reads_R2.fastq.gz \       # 反向reads
    --careful \                  # 仔细模式，减少错误但更慢
    -k 21,33,55,77 \             # k-mer值列表
    -t 16 \                      # CPU线程数
    -m 64 \                      # 最大内存使用量(GB)
    -o spades_output              # 输出目录

# === 使用Flye进行长读段组装 ===
# 适合：PacBio HiFi / Nanopore长读段
flye \
    --pacbio-hifi reads.fastq.gz \  # PacBio HiFi数据（或用--nano-raw指定Nanopore）
    --out-dir flye_output \          # 输出目录
    --threads 16 \                   # CPU线程数
    --genome-size 500m               # 预估基因组大小（500Mb）

# === 使用hifiasm进行HiFi组装 ===
# 2024-2025最推荐的HiFi组装工具
hifiasm \
    -o output_prefix \       # 输出文件前缀
    -t 16 \                  # CPU线程数
    --h1 hi-c_R1.fastq.gz \  # Hi-C数据（可选，用于分型）
    --h2 hi-c_R2.fastq.gz \  # Hi-C数据
    reads.hifi.fastq.gz       # HiFi reads输入

# === 使用MaSuRCA进行混合组装 ===
# 结合二代高准确度 + 三代高连续性
# 第1步：创建配置文件
cat > sr_config.txt << 'EOF'
DATA
PE = pe 300 50 reads_R1.fastq.gz reads_R2.fastq.gz    # 二代双端数据，insert size 300±50
NANOPORE = nanopore.fastq.gz                            # Nanopore长读段
END

PARAMETERS
GRAPH_KMER_SIZE = auto    # 自动选择k-mer大小
USE_LINKING_MATES = 0     # 不使用linking mates
MEGA_READS_ONE_PASS = 0   # 双pass模式，更准确
NUM_THREADS = 16           # CPU线程数
JF_SIZE = 200000000        # Jellyfish哈希表大小
END
EOF

# 第2步：生成组装脚本
masurca sr_config.txt      # 根据配置生成assemble.sh

# 第3步：运行组装
bash assemble.sh           # 执行组装流程

# === 使用QUAST评估组装质量 ===
quast.py \
    assembly.fasta \              # 待评估的组装结果
    -r reference.fasta \          # 参考基因组（可选）
    -o quast_output \             # 输出目录
    -t 8 \                        # CPU线程数
    --large                        # 大基因组模式

# === 关键指标解读 ===
# N50: 越大越好，说明组装连续性好
# L50: 需要多少条序列才能覆盖50%基因组，越小越好
# Total length: 应接近预估基因组大小
# # contigs: 越少越好
# Misassemblies: 错误组装数，越少越好

# === 使用BUSCO评估基因组完整性 ===
busco \
    -i assembly.fasta \                    # 输入组装结果
    -o busco_output \                      # 输出目录
    -m genome \                            # 模式：genome/proteins/transcriptome
    -l eukaryota_odb10 \                   # 谱系数据库（根据物种选择）
    -c 8                                    # CPU线程数

# BUSCO结果解读：
# Complete: >95% 说明组装很完整
# Duplicated: 高比例可能是杂合或多倍体
# Fragmented: 片段化的基因，组装不完整
# Missing: 缺失的基因

# === Pilon纠错（用二代数据纠正三代组装） ===
# 第1步：将二代reads比对回组装结果
bwa index assembly.fasta                                    # 建立索引
bwa mem -t 16 assembly.fasta R1.fq.gz R2.fq.gz | \        # 比对
    samtools sort -@ 8 -o aligned.bam                       # 排序
samtools index aligned.bam                                  # 建立BAM索引

# 第2步：运行Pilon纠错
java -Xmx64G -jar pilon.jar \
    --genome assembly.fasta \    # 原始组装结果
    --frags aligned.bam \        # 比对结果
    --output polished \          # 输出前缀
    --threads 16                  # CPU线程数

# === Purge_dups去除冗余序列 ===
# 适用于杂合基因组，去除重复的等位基因区域
minimap2 -xmap-pb assembly.fasta reads.fastq.gz | \   # 将reads比对到组装
    gzip > aligned.paf.gz                                # 压缩PAF文件

pbcstat aligned.paf.gz           # 统计覆盖度
calcuts PB.stat > cutoffs.txt    # 计算覆盖度阈值

purge_dups \
    -2 \                         # 二次过滤
    -T cutoffs.txt \             # 覆盖度阈值
    -c PB.base.cov \             # 覆盖度文件
    assembly.fasta                # 输入组装结果