生物信息学入门实验项目

第一章 项目介绍与学习目标

1.1 项目背景

二十一世纪以来，高通量测序技术（Next-Generation Sequencing, NGS）的飞速发展与测序成本的指数级下降，使生物学研究全面迈入"大数据"时代。如今，一台桌面级测序仪（如 Illumina MiSeq）即可在数小时内产出上千万条、总量达数亿碱基的 DNA 序列数据。面对如此规模的数据，传统的"手工比对、肉眼检查"早已无能为力，生物信息学（Bioinformatics） 应运而生，成为连接测序数据与生物学发现之间不可或缺的桥梁。

生物信息学是一门交叉学科，它综合运用计算机科学、数学与统计学的方法，对生物数据进行存储、处理、分析与解读。对一名生命科学或数据科学方向的学生而言，掌握生物信息学的核心分析流程，已从"加分项"逐渐演变为"基本功"。

然而，入门生物信息学常常面临两大困难：其一，概念抽象、工具繁多，初学者容易在命令行的汪洋中迷失方向；其二，缺乏一条完整、可复现、数据量适中的主线来串联零散的知识点。本实验项目正是为解决这一痛点而设计。

1.1.1 为什么选择"大肠杆菌重测序"作为入门项目

本项目以大肠杆菌（Escherichia coli）全基因组重测序与变异检测为主线。选择这一主题，是基于以下考量：

• 基因组小、计算友好：大肠杆菌 K-12 MG1655 的参考基因组仅约 4.6 Mb（兆碱基），不到人类基因组的千分之一点五。所有分析步骤（质控、比对、变异检测）都可在普通个人电脑上数分钟至数十分钟内完成，适合课堂环境与入门练习。
• 数据公开、真实可复现：本项目所用的测序数据与参考基因组均来自国际公共数据库（ENA、NCBI），任何学生都可免费下载、复现全部结果。
• 流程完整、代表性强：从原始测序数据到变异位点，覆盖了生物信息学最核心的一条"黄金流程"——数据获取 → 质量控制 → 序列比对 → 比对后处理 → 变异检测 → 变异注释。这条流程的每个环节都是通用技能，可直接迁移到人类、动植物、微生物等其他物种的分析中。
• 结果可解读、有科学内涵：通过对比"同源对照样本"与"野生分离株样本"，学生既能掌握变异检测的技术细节，又能深入理解测序深度、质量控制与假阳性之间的本质关系。

1.2 项目目标

1.2.1 知识目标

1. 理解高通量测序（Illumina）的基本原理与数据特征；
2. 掌握 FASTA、FASTQ、SAM/BAM、VCF 等核心生物信息学文件格式的结构与含义；
3. 理解序列比对（BWA-MEM）与变异检测（bcftools）背后的基本算法思想；
4. 理解 Phred 质量分数、测序深度、覆盖度、假阳性等关键概念；
5. 掌握变异注释（SnpEff）的意义与变异类型的分类方法。

1.2.2 技能目标

1. 熟练使用 Linux 命令行完成文件的查看、压缩、统计等基本操作；
2. 使用 conda/mamba 搭建可复现的生物信息学分析环境；
3. 从公共数据库下载参考基因组与测序数据，并校验数据完整性；
4. 使用 FastQC 与 fastp 进行测序数据质量评估与过滤；
5. 使用 BWA-MEM 将测序读段比对到参考基因组，并使用 samtools 完成排序、去重、统计；
6. 使用 bcftools 进行变异位点检测与过滤，并解读 VCF 结果；
7. 使用 SnpEff 对变异进行功能注释与分类统计；
8. 撰写规范的分析报告，以图表与文字呈现分析结果。

1.2.3 能力目标（可迁移的通用能力）

• 工程化思维：学会通过目录组织、脚本化、记录命令等方式管理一个可复现的分析项目；
• 排错能力：学会阅读错误信息、查阅软件文档与社区资源独立解决问题；
• 数据素养：学会对分析结果保持批判性审视，区分"技术噪音"与"真实信号"。

1.3 项目总览

1.3.1 分析流程

本项目的主干流程如图 1.1 所示，共分为七个模块。整个流程是一个典型的"变异检测流水线"（Variant Calling Pipeline）。

图 1.1：本项目主干分析流程

1.3.2 数据集方案

本项目提供两套真实测序数据集，形成"主实验 + 延伸实验"的递进设计：

表 1.1：项目两套数据集的对比

1.4 预备知识要求

本项目面向本科二、三年级以上或具有同等基础的初学者，建议具备以下先修知识：

• 生物学基础：了解 DNA 双螺旋结构、碱基配对（A-T、C-G）、中心法则（DNA → RNA → 蛋白质）的基本概念；
• 计算机基础：会使用 Linux 基本命令（cd、ls、mkdir、cp、mv），了解绝对路径与相对路径的概念；
• 英语基础：能借助词典阅读软件帮助文档与数据库页面（本项目的命令与报错信息均为英文）。

第二章 实验环境与软件清单

2.1 硬件与系统要求

2.2 所需软件清单

所有软件均通过 conda（或更快的 mamba）统一安装。

2.3 安装与环境配置

2.3.1 安装 Miniconda

# 下载 Miniconda 安装脚本
wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
# 运行安装脚本
bash Miniconda3-latest-Linux-x86_64.sh
# 按提示操作：回车阅读协议 → 输入 yes 接受 → 回车确认安装路径 → 输入 yes 初始化

安装完成后关闭并重新打开终端（或执行 source ~/.bashrc），验证：

conda --version
# 输出形如：conda 24.11.0

2.3.2 配置 conda 源（可选但推荐）

# 添加清华镜像
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda/
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main/
conda config --set show_channel_urls yes

2.3.3 创建独立的分析环境

conda create -n bioinfo -y -c conda-forge -c bioconda \
    mamba fastqc fastp bwa samtools bcftools snpeff tabix

conda activate bioinfo

bwa 2>&1 | head -1
samtools --version | head -1
bcftools --version | head -1

2.3.4 下载 SnpEff 数据库

# 查看可用的大肠杆菌数据库
snpEff databases | grep -i coli | head

# 下载 MG1655 (ASM584v2) 数据库
snpEff download -v ASM584v2

2.4 项目目录结构

mkdir -p ecoli_reseq/{data,ref,results/{qc,align,variants,annotation},scripts,report}
cd ecoli_reseq
tree -L 2   # 若未安装 tree 可省略

第三章 测序技术基础

3.1 DNA 测序技术简史

1977 年，Frederick Sanger 发明了链终止法测序（Sanger sequencing），奠定了现代 DNA 测序的基础。其核心原理是"双脱氧核苷酸终止"。Sanger 测序读长可达 800-1000 bp，准确率极高（>99.9%），但通量低。

2005 年前后，以 Roche 454、Illumina（Solexa）、ABI SOLiD 为代表的第二代测序（NGS）技术相继问世。它们的共同特征是大规模并行化。其中 Illumina 平台凭借高通量、低成本、低错误率的综合优势，成为目前应用最广泛的测序平台。

以 PacBio（SMRT）和 Oxford Nanopore（纳米孔测序）为代表的第三代测序，无需 PCR 扩增即可对单条 DNA 分子进行实时测序，读长可达数十 kb 甚至 Mb 级，但单碱基错误率相对较高。

3.2 Illumina 测序原理

本项目的数据由 Illumina MiSeq 平台产生，其测序过程可概括为"桥式扩增 + 边合成边测序"两大步骤。

3.2.1 文库制备

首先将基因组 DNA 随机打断为约几百 bp 的小片段（fragment），在片段两端连接上已知序列的接头（adapter）。

3.2.2 桥式 PCR 扩增（簇生成）

将文库加载到 flow cell 上，DNA 片段通过接头与芯片表面的互补寡核苷酸结合。随后进行"桥式 PCR"，最终形成由同一模板扩增而来的DNA 簇（cluster）。

3.2.3 边合成边测序（SBS）

1. 与模板互补的 dNTP 被 DNA 聚合酶掺入；
2. 激发荧光，记录每个簇发出的颜色，据此判断该位置掺入的碱基；
3. 切除荧光基团与终止子，进入下一轮合成。

如此循环，每次读出一个碱基，最终获得每个簇的序列，即一条测序读段（read）。

3.2.4 单端测序与双端测序

• 单端测序（Single-end, SE）：只从 DNA 片段的一端开始测序，每条片段产生 1 条 read；
• 双端测序（Paired-end, PE）：先从片段一端测序，再翻转到另一端测序，每条片段产生 2 条 read（R1/R2）。

3.3 测序质量与 Phred 分数

测序仪在读碱基时会根据荧光信号强度等指标，为每个碱基赋予一个质量分数，用于衡量"该碱基被读错"的概率。质量分数通常用 Phred 分数（记作 Q）表示：

Q = −10 log₁₀ P

第四章 核心生物信息学文件格式

4.1 FASTA 格式：序列的存储

FASTA 是最简单的序列存储格式。每条序列由两行或多行组成：第一行以 > 开头，为描述行（header）；后续行为序列本身。

>NC_000913.3 Escherichia coli str. K-12 substr. MG1655, complete genome
AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAGC
TTCTGAACTGGTTACCTGCCGTGAGTAAATTAAAATTTTATTGACTTAGGTCACTAAATACTTTAACCAA

4.2 FASTQ 格式：测序读段与质量

FASTQ 是存储测序读段及其质量分数的标准格式，每条 read 占据四行：

1. 第 1 行（标识行）：以 @ 开头，包含 read 的唯一标识符；
2. 第 2 行（序列行）：碱基序列（A/T/C/G/N）；
3. 第 3 行（分隔行）：以 + 开头；
4. 第 4 行（质量行）：与第 2 行等长的质量字符。

@ERR15404827.1 1/1
AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATT
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

4.3 SAM/BAM 格式：比对结果的存储

SAM（Sequence Alignment/Map）是存储"读段比对到参考基因组"结果的文本格式；BAM 是其二进制压缩版本。

4.3.1 SAM 的 11 个必选字段

4.3.2 CIGAR 常用操作符

4.4 VCF/BCF 格式：变异位点的存储

VCF（Variant Call Format）是存储遗传变异（SNP、插入缺失等）的标准格式。前 8 个固定字段含义：

FORMAT 与样本列常用字段：

• GT：基因型（0/1 杂合，1/1 纯合变异）；
• DP：该位点的测序深度；
• AD：各等位基因的 read 深度。

第五章 核心算法原理

5.1 序列比对算法

5.1.1 BWT 与 FM-index

BWA 的核心是基于 Burrows-Wheeler 变换（BWT） 和 FM-index 构建的全文索引。其思想：

1. 将参考基因组的所有"循环移位"按字典序排序，取每行最后一列，得到 BWT 字符串。这个变换是可逆的；
2. FM-index 在 BWT 上支持"反向搜索（backward search）"，可在 O(m) 时间内（m 为 read 长度）快速定位某个短序列在参考基因组中的所有出现位置；
3. 比对一条 read 时，BWA 从 read 末端开始逐碱基"反向搜索"，一旦发现不匹配，立即回溯。

5.1.2 seed-and-extend 策略

BWA-MEM 采用"种子-扩展"（seed-and-extend）策略：

1. 种子查找：从 read 中提取若干短的"种子"，利用 FM-index 快速找到候选位置；
2. 链化（chaining）：将同一条 read 的多个种子候选位置"串联"成可能的一致比对；
3. 扩展（extension）：对每个候选位置，用 Smith-Waterman 局部比对算法向两端扩展。

5.1.3 比对质量 MAPQ

MAPQ（mapping quality）是 BWA 为每条 read 赋予的"比对可信度"（Phred 尺度）：

MAPQ = −10 log₁₀ P(比对位置错误)

当一条 read 能比对到多个位置（如重复区）时，MAPQ 会很低（如 0）。

5.2 变异检测算法

本项目使用的 bcftools mpileup + call 采用"堆积 → 基因型似然 → 变异判定"三步策略。

5.2.1 pileup 堆积

mpileup 将 BAM 文件中覆盖同一基因组位置的所有 read 的碱基"堆积"在一起，记录每个位置 A/C/G/T 各有多少条 read 支持。

5.2.2 计算基因型似然

对每个候选位点，基于观测到的碱基堆积，计算不同基因型假设下的似然（likelihood）。

5.2.3 贝叶斯后验与输出

实际算法会结合先验概率（如 SNP 的先验概率约 10⁻³）与似然，用贝叶斯公式计算后验概率：

QUAL = −10 log₁₀ P(该位置不是变异)

5.3 测序深度与覆盖度

测序深度 = 测序数据总量 / 基因组大小。

第六章 模块 1：数据获取

6.1 背景知识：公共生物信息学数据库

全球三大核酸序列数据库——GenBank（NCBI）、ENA（EBI）、DDBJ——共同组成国际核苷酸序列数据库协作联盟（INSDC）：

6.2 实验步骤

6.2.1 步骤 1：创建目录结构

mkdir -p ~/ecoli_reseq/{data,ref,results/{qc,align,variants,annotation},scripts,report}
cd ~/ecoli_reseq

6.2.2 步骤 2：下载参考基因组

cd ref

wget -c https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.fna.gz
wget -c https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.gff.gz

gunzip GCF_000005845.2_ASM584v2_genomic.fna.gz

ls -lh

6.2.3 步骤 3：查看参考基因组

grep "^>" GCF_000005845.2_ASM584v2_genomic.fna
grep -v "^>" GCF_000005845.2_ASM584v2_genomic.fna | tr -d '\n' | wc -c

预期输出：

>NC_000913.3 Escherichia coli str. K-12 substr. MG1655, complete genome
4641652

6.2.4 步骤 4：下载测序数据

cd ../data

wget -c https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR154/027/ERR15404827/ERR15404827_1.fastq.gz
wget -c https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR154/027/ERR15404827/ERR15404827_2.fastq.gz

ls -lh

6.2.5 步骤 5：数据完整性校验

# 计算本地文件的 MD5 值
md5sum ERR15404827_1.fastq.gz ERR15404827_2.fastq.gz

# 获取 ENA 官方 MD5 校验值
curl -s "https://www.ebi.ac.uk/ena/portal/api/filereport?accession=ERR15404827&result=read_run&fields=run_accession,fastq_md5&format=tsv"

将本地计算的 MD5 与 ENA 官方提供的 MD5 逐字符比对，完全一致即表示文件完整无损。

6.3 思考题

第七章 模块 2：测序数据质量控制

7.1 背景知识：为什么要做质量控制？

原始测序数据中存在的"噪音"包括：测序错误、接头污染、低质量碱基、PCR 重复。

7.2 实验步骤

7.2.1 步骤 1：使用 FastQC 评估

fastqc -t 4 -o results/qc/ \
  data/ERR15404827_1.fastq.gz \
  data/ERR15404827_2.fastq.gz

生成 .html 报告与 .zip 原始数据。

7.2.2 步骤 2：解读 FastQC 报告

重点查看：碱基质量、碱基 GC 含量、接头污染、重复水平、重复序列等指标。FastQC 用绿/黄/红三色标记各项指标。

7.2.3 步骤 3：使用 fastp 过滤与修剪

fastp \
  -i data/ERR15404827_1.fastq.gz \
  -I data/ERR15404827_2.fastq.gz \
  -o results/qc/ERR15404827_1.clean.fastq.gz \
  -O results/qc/ERR15404827_2.clean.fastq.gz \
  --html results/qc/fastp.html \
  --json results/qc/fastp.json \
  --thread 4 \
  --qualified_quality_phred 20 \
  --unqualified_percent_limit 40 \
  --length_required 50 \
  --cut_front --cut_tail

7.2.4 步骤 4：查看 fastp 报告与统计

python3 -c "
import json
d = json.load(open('results/qc/fastp.json'))
s = d['summary']
for k in ['before_filtering','after_filtering']:
    x = s[k]
    print(k, '-> total_reads:', x['total_reads'],
          '| q30_rate: %.2f%%' % (x['q30_rate']*100))
"

7.3 思考题

第八章 模块 3：序列比对

8.1 背景知识：什么是序列比对

序列比对是生物信息学最核心的操作之一。任务是：给定一条短测序读段和一个参考基因组，确定该读段最可能来自参考基因组的哪个位置。

8.2 实验步骤

8.2.1 步骤 1：为参考基因组建立索引

bwa index ref/GCF_000005845.2_ASM584v2_genomic.fna
samtools faidx ref/GCF_000005845.2_ASM584v2_genomic.fna

8.2.2 步骤 2：使用 BWA-MEM 进行比对

bwa mem -t 4 \
  ref/GCF_000005845.2_ASM584v2_genomic.fna \
  results/qc/ERR15404827_1.clean.fastq.gz \
  results/qc/ERR15404827_2.clean.fastq.gz \
  > results/align/ERR15404827.sam

8.2.3 步骤 3：查看 SAM 文件

samtools view -H results/align/ERR15404827.sam | head -10
samtools view results/align/ERR15404827.sam | head -5

8.2.4 步骤 4：评估比对结果

samtools flagstat results/align/ERR15404827.sam

关注比对率（mapped）与正确配对比率（properly paired）。对同源重测序，比对率通常应 >99%，properly paired >95%。

8.3 思考题

第九章 模块 4：比对后处理

9.1 背景知识

BWA 输出的 SAM 存在三个问题：文件巨大且低效、无序、含 PCR 重复。

9.2 实验步骤

9.2.1 步骤 1：SAM 转换为 BAM

samtools view -bS results/align/ERR15404827.sam \
  -o results/align/ERR15404827.bam

9.2.2 步骤 2：按坐标排序

samtools sort -@ 4 results/align/ERR15404827.bam \
  -o results/align/ERR15404827.sorted.bam

9.2.3 步骤 3：标记 mate 信息

samtools fixmate -m results/align/ERR15404827.sorted.bam \
  results/align/ERR15404827.fixmate.bam

9.2.4 步骤 4：再次排序

samtools sort -@ 4 results/align/ERR15404827.fixmate.bam \
  -o results/align/ERR15404827.fixmate.sorted.bam

9.2.5 步骤 5：去除 PCR 重复

samtools markdup results/align/ERR15404827.fixmate.sorted.bam \
  results/align/ERR15404827.dedup.bam

9.2.6 步骤 6：建立 BAM 索引

samtools index results/align/ERR15404827.dedup.bam

9.2.7 步骤 7：比对统计

samtools flagstat results/align/ERR15404827.dedup.bam
samtools stats results/align/ERR15404827.dedup.bam > results/align/stats.txt
samtools coverage results/align/ERR15404827.dedup.bam

9.2.8 步骤 8：计算深度分布

samtools depth -a results/align/ERR15404827.dedup.bam \
  > results/align/depth.txt

awk '{s+=$3; if($3>0) c++} END{print "平均深度(含0):", s/NR; print "覆盖位点平均深度:", s/c; print "覆盖度:", c/NR*100"%"}' \
  results/align/depth.txt

9.3 思考题

第十章 模块 5：变异检测

10.1 背景知识：变异检测的两步流程

本项目使用的 bcftools 将变异检测拆分为两个可组合的步骤：

1. mpileup（堆积）：将覆盖同一位置的 read 堆积，统计 A/C/G/T 支持数与质量；
2. call（判定）：计算不同基因型假设的似然，结合先验概率做贝叶斯判定。

10.2 实验步骤

10.2.1 步骤 1：变异检测

bcftools mpileup --ploidy 1 \
  -f ref/GCF_000005845.2_ASM584v2_genomic.fna \
  results/align/ERR15404827.dedup.bam | \
bcftools call -mv -Ov \
  -o results/variants/ERR15404827.raw.vcf

10.2.2 步骤 2：查看原始 VCF

grep -vc "^#" results/variants/ERR15404827.raw.vcf
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%QUAL\t%DP\n' \
  results/variants/ERR15404827.raw.vcf | head -10

10.2.3 步骤 3：变异质量过滤

bcftools filter \
  -i 'QUAL >= 30 && DP >= 10' \
  results/variants/ERR15404827.raw.vcf \
  -o results/variants/ERR15404827.filtered.vcf

grep -vc "^#" results/variants/ERR15404827.filtered.vcf

10.2.4 步骤 4：变异统计

bcftools stats results/variants/ERR15404827.filtered.vcf \
  > results/variants/stats.txt
grep -E "SN|number of" results/variants/stats.txt | head -30

10.2.5 步骤 5：解读主实验的"近乎零变异"结果

10.3 思考题

第十一章 模块 6：变异注释与功能解读

11.1 背景知识：为什么要注释变异

变异检测得到的 VCF 只是一个"位置 + 碱基变化"的列表。变异注释就是为每个变异补充其功能信息。

11.2 实验步骤

11.2.1 步骤 1：运行 SnpEff 注释

snpEff -v ASM584v2 \
  results/variants/ERR15404827.filtered.vcf \
  > results/annotation/ERR15404827.ann.vcf

11.2.2 步骤 2：查看注释结果

bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%INFO/ANN\n' \
  results/annotation/ERR15404827.ann.vcf | head -10

11.2.3 步骤 3：理解 ANN 字段

11.2.4 步骤 4：变异类型分类统计

bcftools query -f '%INFO/ANN\n' \
  results/annotation/ERR15404827.ann.vcf \
  | cut -d'|' -f2 | sort | uniq -c | sort -rn

11.3 思考题

第十二章 模块 7：结果汇总与延伸实验

12.1 步骤 1：汇总主实验关键结果

12.2 步骤 2：撰写分析报告

分析报告结构：标题与摘要、数据与方法、结果、讨论、结论。

12.3 步骤 3：延伸实验——野生分离株分析

12.3.1 延伸实验数据下载

wget -c https://ftp.sra.ebi.ac.uk/vol1/fastq/DRR063/DRR063436/DRR063436_1.fastq.gz
wget -c https://ftp.sra.ebi.ac.uk/vol1/fastq/DRR063/DRR063436/DRR063436_2.fastq.gz

md5sum DRR063436_1.fastq.gz DRR063436_2.fastq.gz

12.3.2 延伸实验完整流程

流程与主实验完全一致，只需将文件名 ERR15404827 替换为 DRR063436。

fastp \
  -i data/DRR063436_1.fastq.gz \
  -I data/DRR063436_2.fastq.gz \
  -o results/qc/DRR063436_1.clean.fastq.gz \
  -O results/qc/DRR063436_2.clean.fastq.gz \
  --html results/qc/fastp_DRR063436.html \
  --json results/qc/fastp_DRR063436.json \
  --thread 4 \
  --qualified_quality_phred 20 \
  --unqualified_percent_limit 40 \
  --length_required 40

12.4 步骤 4：对比分析与讨论

12.5 项目总结

附录 A 软件与环境清单汇总

A.1 完整软件清单

A.2 一键式环境配置

A.2.1 方式一：逐条命令安装

conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda/
conda config --set show_channel_urls yes

conda create -n bioinfo -y -c conda-forge -c bioconda \
    mamba fastqc fastp bwa samtools bcftools snpeff tabix

conda activate bioinfo
snpEff download -v ASM584v2

A.2.2 方式二：使用环境文件

将以下内容保存为 environment.yml：

name: bioinfo
channels:
  - conda-forge
  - bioconda
dependencies:
  - mamba
  - fastqc=0.12.1
  - fastp=0.23.4
  - bwa=0.7.18
  - samtools=1.20
  - bcftools=1.20
  - snpeff=5.2
  - tabix=1.20

conda env create -f environment.yml
conda activate bioinfo
snpEff download -v ASM584v2

A.3 Windows 用户：安装 WSL2

Windows 10/11 用户可通过 WSL2 获得完整的 Linux 环境：

# 安装 WSL 与 Ubuntu（在 Windows PowerShell 中以管理员身份运行）
wsl --install -d Ubuntu-22.04
# 重启后，Ubuntu 会自动启动，设置用户名与密码

附录 B 常见问题与排错

附录 C 术语表

附录 D 全部命令速查

D.1 模块 1：数据获取

cd ref
wget -c https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.fna.gz
wget -c https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.gff.gz
gunzip GCF_000005845.2_ASM584v2_genomic.fna.gz

cd ../data
wget -c https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR154/027/ERR15404827/ERR15404827_1.fastq.gz
wget -c https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR154/027/ERR15404827/ERR15404827_2.fastq.gz

md5sum *.fastq.gz

D.2 模块 2：质量控制

fastqc -t 4 -o results/qc/ data/*.fastq.gz
fastp -i data/ERR15404827_1.fastq.gz -I data/ERR15404827_2.fastq.gz \
  -o results/qc/ERR15404827_1.clean.fastq.gz -O results/qc/ERR15404827_2.clean.fastq.gz \
  --html results/qc/fastp.html --json results/qc/fastp.json \
  --thread 4 --qualified_quality_phred 20 --unqualified_percent_limit 40 \
  --length_required 50 --cut_front --cut_tail

D.3 模块 3：序列比对

bwa index ref/GCF_000005845.2_ASM584v2_genomic.fna
samtools faidx ref/GCF_000005845.2_ASM584v2_genomic.fna
bwa mem -t 4 ref/GCF_000005845.2_ASM584v2_genomic.fna \
  results/qc/ERR15404827_1.clean.fastq.gz results/qc/ERR15404827_2.clean.fastq.gz \
  > results/align/ERR15404827.sam

D.4 模块 4：比对后处理

samtools view -bS results/align/ERR15404827.sam -o results/align/ERR15404827.bam
samtools sort -@ 4 results/align/ERR15404827.bam -o results/align/ERR15404827.sorted.bam
samtools fixmate -m results/align/ERR15404827.sorted.bam results/align/ERR15404827.fixmate.bam
samtools sort -@ 4 results/align/ERR15404827.fixmate.bam -o results/align/ERR15404827.fixmate.sorted.bam
samtools markdup results/align/ERR15404827.fixmate.sorted.bam results/align/ERR15404827.dedup.bam
samtools index results/align/ERR15404827.dedup.bam
samtools flagstat results/align/ERR15404827.dedup.bam
samtools coverage results/align/ERR15404827.dedup.bam

D.5 模块 5：变异检测

bcftools mpileup --ploidy 1 -f ref/GCF_000005845.2_ASM584v2_genomic.fna \
  results/align/ERR15404827.dedup.bam | \
bcftools call -mv -Ov -o results/variants/ERR15404827.raw.vcf
bcftools filter -i 'QUAL >= 30 && DP >= 10' \
  results/variants/ERR15404827.raw.vcf -o results/variants/ERR15404827.filtered.vcf
bcftools stats results/variants/ERR15404827.filtered.vcf > results/variants/stats.txt

D.6 模块 6：变异注释

snpEff -v ASM584v2 results/variants/ERR15404827.filtered.vcf > results/annotation/ERR15404827.ann.vcf
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%INFO/ANN\n' results/annotation/ERR15404827.ann.vcf

附录 E 思考题参考答案

E.1 模块 1

1. 为什么用 .gz 压缩？ FASTQ/FASTA 文件体积庞大，gzip 压缩可节省 70-80% 存储与传输成本；且 gzip 是流式压缩，工具可直接读取 .gz 文件。
2. R1/R2 的含义 双端测序从同一 DNA 片段的两端各测一条 read，分别存入 R1、R2 文件。若遗漏 R2，则退化为单端数据，比对准确性下降。
3. MD5 不一致的可能原因 文件传输损坏、下载不完整、文件被篡改、或下载了错误版本的文件。

E.2 模块 2

1. 说明测序循环后期信号衰减、错误率升高。fastp 的 --cut_tail 会从末端滑动窗口修剪低质量碱基。
2. 过滤会丢弃低质量/过短的 read，故总数减少。损失 >30% 说明数据整体质量差、或过滤参数过严。
3. 当插入片段短于测序读长时，测序"读穿"了整个片段，继续读到另一端的接头序列。

E.3 模块 3

1. 为什么先建索引？ 索引把"逐位置扫描比对"转化为"常数级快速查询"，使得数百万条 read 比对成为可能。
2. FLAG=99 = 1+2+32+64，即：双端测序、与 mate 均正确比对、mate 在负链、是 read1。
3. MAPQ=0 表示比对位置不唯一（可能比对到多个位置，如重复区）。

E.4 模块 4

1. 按坐标排序后，同一位置的 read 在文件中连续存放，变异检测可以顺序"堆积"同一位置的 read，避免随机访问带来的巨大开销。
2. PCR 重复是同一 DNA 分子的多个扩增拷贝，信息冗余且会放大测序错误；生物学重复是独立来源的样本。变异检测需独立证据，PCR 重复只会放大错误。
3. 用"数据量/基因组大小"估算为约 124 Mb / 4.64 Mb ≈ 27×。实际值略低的可能原因：部分 read 被过滤/去重、部分 read 比对到参考基因组之外的序列。

E.5 模块 5

1. 二者都是 Phred 尺度的质量分数，但对象不同：FASTQ 的碱基质量衡量"该碱基被测序仪读错"的概率；VCF 的 QUAL 衡量"该位点不是变异（判定错误）"的概率。
2. QUAL 反映变异判定的统计置信度，DP 反映证据量。只设 QUAL 会放过"高置信但证据薄弱"的位点。
3. 更严的过滤会进一步减少变异数，去掉更多疑似假阳性，但也可能漏掉真实变异（尤其低深度区域）。

E.6 模块 6

1. 遗传密码具有简并性（多个密码子编码同一氨基酸），故密码子第三位常发生同义替换而不改变氨基酸。
2. 无义突变影响通常更大：它导致蛋白提前截断，几乎必然破坏功能。
3. 基因间区不编码蛋白，变异通常不直接影响蛋白功能，因而更"容忍"、更易保留。

附录 F 数据集信息卡与报告模板

F.1 数据集信息卡

F.1.1 主实验数据集：ERR15404827

F.1.2 延伸实验数据集：DRR063436

F.1.3 参考基因组

F.2 分析报告模板

# 大肠杆菌重测序数据分析报告

## 1. 摘要
（一段话概括：目的、数据、方法、主要结论）

## 2. 数据与方法
- 数据集：ERR15404827（主）/ DRR063436（延伸）
- 参考基因组：GCF_000005845.2 (NC_000913.3)
- 软件版本：FastQC 0.12.1、fastp 0.23.4、BWA 0.7.18、
  samtools 1.20、bcftools 1.20、SnpEff 5.2
- 关键参数：质量过滤 Q>=20；比对 BWA-MEM 默认；
  变异检测 ploidy=1；过滤 QUAL>=30 && DP>=10

## 3. 结果
- 3.1 质控统计
- 3.2 比对统计
- 3.3 变异检测结果
- 3.4 变异注释结果

## 4. 讨论
- 同源对照接近零变异的科学意义
- 深度对变异检测的影响（对比主/延伸实验）
- 局限性与改进方向

## 5. 结论
（概括最重要的发现与收获）