609_宏基因组组装评估QUAST¶

# 创建新环境（避免依赖冲突）
conda create -n quast_env python=3.8  # 创建Python3.8环境
conda activate quast_env  # 激活环境

# 安装QUAST（包含MetaQUAST）
conda install -c bioconda quast  # 从bioconda频道安装

# 验证安装
quast --version  # 查看版本号，当前最新5.3.0
metaquast.py --help  # 查看MetaQUAST帮助信息

pip install quast  # 通过pip安装

git clone https://github.com/ablab/quast.git  # 克隆仓库
cd quast  # 进入目录
python setup.py install  # 安装

# 只评估组装本身的统计信息
quast contigs.fasta \
    -o quast_output \          # 输出目录
    -t 8                       # 使用8个线程

# 有参考基因组时，能计算错配率、基因覆盖率等
quast contigs.fasta \
    -r reference.fasta \       # 参考基因组
    -g genes.gff \             # 参考基因组的基因注释（可选）
    -o quast_with_ref \        # 输出目录
    -t 8 \                     # 线程数
    --min-contig 500           # 最短contig长度过滤（默认500bp）

# 同时比较不同组装器的结果
quast megahit_contigs.fasta \
    spades_contigs.fasta \
    idba_contigs.fasta \
    -r reference.fasta \       # 参考基因组
    -l "MEGAHIT,SPAdes,IDBA" \ # 给每个组装起标签
    -o compare_assemblers \    # 输出目录
    -t 16                      # 线程数

# 提供已知参考基因组
metaquast.py metagenome_contigs.fasta \
    -R ref_genome1.fasta,ref_genome2.fasta \  # 多个参考基因组
    -o metaquast_output \      # 输出目录
    -t 16 \                    # 线程数
    --max-ref-number 50        # 最多使用50个参考基因组

# MetaQUAST自动用BLAST搜索SILVA数据库找参考
metaquast.py metagenome_contigs.fasta \
    -o metaquast_auto \        # 输出目录
    -t 16 \                    # 线程数
    --max-ref-number 50        # 最多下载50个参考基因组

# 比较MEGAHIT和metaSPAdes的组装结果
metaquast.py megahit_contigs.fa spades_contigs.fa \
    -l "MEGAHIT,metaSPAdes" \  # 标签
    -R references/ \           # 参考基因组目录
    -o meta_compare \          # 输出目录
    -t 16 \                    # 线程数
    --min-identity 90          # 最低比对一致性（默认90%）

# 查看文本报告
cat quast_output/report.txt

# 用浏览器打开HTML报告
firefox quast_output/report.html  # Linux
# 或
open quast_output/icarus.html     # macOS

import pandas as pd  # 导入pandas

# 读取QUAST的TSV报告
df = pd.read_csv("quast_output/report.tsv", sep="\t")  # 读取报告
print(df.head(20))  # 打印前20行

# 提取关键指标
n50 = df[df["Assembly"] == "N50"].iloc[0, 1]  # 获取N50值
total_len = df[df["Assembly"] == "Total length"].iloc[0, 1]  # 获取总长度
print(f"N50: {n50}")  # 打印N50
print(f"Total length: {total_len}")  # 打印总长度

#!/bin/bash
# 完整的宏基因组组装评估流程

# === 第1步：组装（用MEGAHIT和metaSPAdes分别组装）===
megahit -1 reads_R1.fq.gz -2 reads_R2.fq.gz \
    -o megahit_out \           # MEGAHIT输出目录
    -t 16                      # 16线程

metaspades.py -1 reads_R1.fq.gz -2 reads_R2.fq.gz \
    -o spades_out \            # SPAdes输出目录
    -t 16                      # 16线程

# === 第2步：用MetaQUAST评估 ===
metaquast.py \
    megahit_out/final.contigs.fa \     # MEGAHIT组装结果
    spades_out/contigs.fasta \         # SPAdes组装结果
    -l "MEGAHIT,metaSPAdes" \          # 标签
    -o assembly_evaluation \           # 输出目录
    -t 16 \                            # 线程数
    --min-contig 1000 \                # 只评估>1kb的contig
    --gene-finding                     # 启用基因预测

# === 第3步：查看结果 ===
echo "=== 组装评估报告 ==="
cat assembly_evaluation/report.txt     # 打印文本报告

# QUAST常用参数
--min-contig 500       # 最短contig过滤长度（默认500）
--threads / -t         # 线程数
--gene-finding         # 启用基因预测（GeneMark/MetaGeneMark）
--large                # 大基因组模式（QUAST-LG）
--fragmented           # 参考基因组不完整时使用
--no-plots             # 不生成图片（节省时间）
--split-scaffolds      # 在N处断开scaffold再评估
--fast                 # 快速模式（跳过部分分析）

# MetaQUAST特有参数
--max-ref-number 50    # 最多参考基因组数
--min-identity 90      # 最低比对一致性（默认90%）
--references-list      # 参考基因组列表文件
--unique-mapping       # 每个contig只比对到最佳参考