基因组组装

基因组组装（Genome Assembly）是生物信息学的一个核心任务，旨在将从高通量测序技术（如Illumina、PacBio、Nanopore等）获得的短序列片段（reads）拼接成完整的基因组序列。这个过程复杂且需要考虑不同的技术、算法和错误校正。

基因组组装流程

基因组组装流程通常包括以下几个步骤：

1. 数据准备与质量控制

• 目标：确保原始测序数据的高质量。

• 方法：

• 去除低质量reads。

• 过滤掉可能的污染序列（如细菌或病毒污染）。

• 去除接头序列（adapters）。

常用工具：FastQC（质量评估）、Trimmomatic（去低质量片段）、BBMap（污染检测）。

2. 重复序列和k-mer分析

• 目标：估计基因组大小、复杂度、重复序列比例。

• 方法：基于k-mer统计分析，计算覆盖深度和重复率。

常用工具：Jellyfish（k-mer计数）、GenomeScope（基因组复杂性估计）。

3. 基因组组装

基因组组装分为de novo组装和参考组装两类。

(1) de novo组装

• 目标：在没有参考基因组的情况下组装完整基因组。

• 方法：

• 重叠图算法（Overlap-Layout-Consensus, OLC）：适用于长读长（PacBio, Nanopore）。

• de Bruijn图算法：适用于短读长（Illumina）。

• 混合算法：结合长短读长的优点。

常用工具：

• OLC方法：Canu, Flye, wtdbg2

• de Bruijn图：SPAdes, Velvet

(2) 参考组装

• 目标：基于已有参考基因组，对序列片段进行比对和扩展。

• 方法：通过序列比对工具，将短reads比对到参考序列并重建基因组。

常用工具：BWA, Bowtie2, SAMtools

4. 组装后处理

• 目标：优化基因组组装结果，使之更接近真实基因组。

• 步骤：

• 工具：GapCloser

• 指标：N50、L50、基因组完整性、错误率。

• 工具：QUAST, BUSCO

• 工具：Pilon, Racon

1. 错误校正：通过对比reads修正组装中的错配或插入缺失。

2. 组装评估：评估基因组组装质量。

3. Gap填充：填补组装中出现的缺口。

5. 注释与下游分析

组装完成后，还需要进行基因组注释（如基因功能预测）、结构变异分析和进化分析等。

基因组组装分析

1. 数据准备与质量控制

import subprocess
# 定义原始数据路径raw_reads = "raw_reads.fastq"clean_reads = "clean_reads.fastq"
def quality_control(input_file, output_file):    """    使用 Trimmomatic 对测序数据进行质量控制    """    cmd = (        f"trimmomatic SE {input_file} {output_file} "        "LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 MINLEN:36"    )    subprocess.run(cmd, shell=True, check=True)    print(f"质量控制完成，输出文件: {output_file}")
# 执行质量控制quality_control(raw_reads, clean_reads)

2. k-mer 分析

def kmer_analysis(input_file, k=31, output_file="kmer_counts.jf"):    """    使用 Jellyfish 进行 k-mer 计数    """    cmd = f"jellyfish count -m {k} -s 100M -C {input_file} -o {output_file}"    subprocess.run(cmd, shell=True, check=True)    print(f"k-mer 分析完成，输出文件: {output_file}")
# 执行 k-mer 分析kmer_analysis(clean_reads)
def genome_complexity_estimation(kmer_counts, output_dir="genome_scope"):    """    使用 GenomeScope 进行基因组复杂性估计    """    cmd = f"jellyfish histo {kmer_counts} | genomescope2 -o {output_dir} -k 31"    subprocess.run(cmd, shell=True, check=True)    print(f"基因组复杂性分析完成，结果保存在: {output_dir}")
# 执行基因组复杂性分析genome_complexity_estimation("kmer_counts.jf")

3. 基因组组装

def genome_assembly(reads, output_dir="assembly_output"):    """    使用 SPAdes 进行 de novo 基因组组装    """    cmd = f"spades.py -s {reads} -o {output_dir}"    subprocess.run(cmd, shell=True, check=True)    print(f"基因组组装完成，结果保存在: {output_dir}")
# 执行基因组组装genome_assembly(clean_reads)

4. 组装后错误校正

def error_correction(assembly_file, reads, corrected_file="corrected_assembly.fasta"):    """    使用 Pilon 对组装结果进行错误校正    """    cmd = (        f"pilon --genome {assembly_file} --frags {reads} --output {corrected_file} "        "--changes --verbose"    )    subprocess.run(cmd, shell=True, check=True)    print(f"错误校正完成，修正后基因组文件: {corrected_file}")
# 执行错误校正error_correction("assembly_output/contigs.fasta", clean_reads)

5. 基因组组装评估

def evaluate_assembly(assembly_file, reference_file=None):    """    使用 QUAST 评估基因组组装质量    """    if reference_file:        cmd = f"quast.py {assembly_file} -r {reference_file} -o quast_results"    else:        cmd = f"quast.py {assembly_file} -o quast_results"    subprocess.run(cmd, shell=True, check=True)    print("基因组组装评估完成，结果保存在: quast_results")
# 执行评估（提供参考基因组时）evaluate_assembly("corrected_assembly.fasta", reference_file="reference_genome.fasta")

6. 可视化结果

以下代码以绘制 N50 值和 contig 长度分布图为例：

import matplotlib.pyplot as plt
def parse_quast_results(quast_report="quast_results/report.txt"):    """    解析 QUAST 生成的报告文件，提取 N50 和 contig 长度分布信息    """    n50 = 0    contig_lengths = []
    with open(quast_report, "r") as f:        for line in f:            if line.startswith("N50"):                n50 = int(line.split()[-1])            if line.startswith(">="):  # 解析 contig 长度分布                contig_lengths.append(int(line.split()[0][2:]))
    return n50, contig_lengths
def plot_results(n50, contig_lengths):    """    绘制 N50 值和 contig 长度分布图    """    plt.figure(figsize=(10, 6))
    # 绘制 N50    plt.subplot(1, 2, 1)    plt.bar(["N50"], [n50], color="skyblue")    plt.ylabel("Length (bp)")    plt.title("N50 Value")
    # 绘制 contig 长度分布    plt.subplot(1, 2, 2)    plt.hist(contig_lengths, bins=30, color="lightgreen")    plt.xlabel("Contig Length (bp)")    plt.ylabel("Frequency")    plt.title("Contig Length Distribution")
    plt.tight_layout()    plt.show()
# 解析并可视化n50, contig_lengths = parse_quast_results()plot_results(n50, contig_lengths)

代码运行注意事项

1. 工具路径配置：

• 确保已安装所需工具（Trimmomatic, Jellyfish, SPAdes, Pilon, QUAST等）。

• 将工具的可执行文件路径加入环境变量或在代码中显式指定。

2. 硬件要求：

• 基因组组装和复杂性分析对内存需求较高，建议使用高性能计算资源。

3. 输入文件格式：

• 确保输入文件格式符合工具要求（如 FASTQ 格式的 reads 和 FASTA 格式的参考基因组）。

4. 数据量：

• 小基因组（如细菌）可以使用以上流程，复杂生物（如人类、植物）的基因组可能需要调整工具和参数。