生信核心数据格式解析手册：FASTQ/BAM/VCF/GFF 处理原理与实操

一、FASTQ 格式：原始测序数据标准格式

1. 基本结构与原理

FASTQ 是存储生物序列及其质量评分的标准格式，每条记录占 4 行：

• 第 1 行：以@开头的序列标识符（含测序运行、簇等信息）
• 第 2 行：核苷酸序列（A、T、C、G、N）
• 第 3 行：以+开头（可重复 ID）
• 第 4 行：质量值（Phred+33/64 编码的 ASCII 字符）

质量值原理：

• Phred 质量值：Q = -10log₁₀(P)，P 为碱基错误概率
• Phred+33 编码：字符 ASCII 值 - 33 = 质量得分（Illumina 1.8 + 标准）
• 常见值：Q30（错误率 0.1%，对应字符 '?'），Q20（错误率 1%，对应 '!'）

2. 处理实操与代码模板

(1) Python：使用 Biopython 解析与质控

python

运行

# 安装：pip install biopython
from Bio import SeqIO
from Bio.SeqIO.QualityIO import FastqGeneralIterator

# 基础读取：统计序列数量
record_count = 0
with open("input.fastq", "r") as handle:
    for record in SeqIO.parse(handle, "fastq"):
        record_count += 1
print(f"总序列数: {record_count}")

# 质量过滤：保留平均质量≥20的reads
min_qual = 20
filtered_records = []
with open("input.fastq", "r") as in_handle, open("output.fastq", "w") as out_handle:
    for title, seq, qual in FastqGeneralIterator(in_handle):
        avg_qual = sum(ord(c) - 33 for c in qual) / len(qual)
        if avg_qual >= min_qual:
            out_handle.write(f"@{title}\n{seq}\n+\n{qual}\n")

(2) 命令行工具：fastp 实现高效质控

bash

运行

# 安装：conda install fastp
fastp -i input.fastq -o output.fastq \
      --cut_mean_quality 20 \   # 平均质量<20的窗口截断
      --length_required 30 \  # 保留长度≥30bp的reads
      --detect_adapter_for_pe  # 自动检测并去除接头

二、BAM 格式：比对结果的高效存储

1. 基本结构与原理

BAM 是 SAM（序列比对 / 映射格式）的二进制压缩形式，具有三大优势：

• 存储效率高：比 SAM 小 70-90%
• 读写速度快：适合大规模数据处理
• 支持索引：可快速定位特定基因组区域

核心组成：

• 文件头：记录参考序列信息
• 比对记录：每条包含 read ID、染色体、起始位置、比对质量 (MAPQ)、CIGAR 字符串（描述比对方式）等

2. 处理实操与代码模板

(1) Python：使用 pysam 操作 BAM 文件

python

运行

# 安装：pip install pysam
import pysam

# 读取BAM文件
bamfile = pysam.AlignmentFile("alignment.bam", "rb")

# 统计比对信息
total_reads = bamfile.mapped + bamfile.unmapped
mapped_rate = bamfile.mapped / total_reads * 100
print(f"总reads: {total_reads}, 比对上: {bamfile.mapped} ({mapped_rate:.2f}%)")

# 提取特定区域比对（chr1:1000-2000）
region_reads = []
for read in bamfile.fetch("chr1", 1000, 2000):
    region_reads.append(read)
print(f"区域chr1:1000-2000包含{len(region_reads)}条reads")

# 保存过滤后的BAM（仅保留MAPQ≥30的reads）
filtered_bam = pysam.AlignmentFile("filtered.bam", "wb", header=bamfile.header)
for read in bamfile:
    if read.mapping_quality >= 30:
        filtered_bam.write(read)

(2) 命令行工具：samtools 实现 BAM 高效处理

bash

运行

# 安装：conda install samtools

# 查看BAM统计信息
samtools stats alignment.bam > stats.txt

# 提取高质量比对（MAPQ≥30）
samtools view -bh -q 30 alignment.bam > high_qual.bam

# 按染色体排序并建立索引
samtools sort alignment.bam -o sorted.bam
samtools index sorted.bam

三、VCF 格式：基因组变异的标准表示

1. 基本结构与原理

VCF（变异调用格式）用于存储基因组变异信息，包括：

• 文件头：以##开头，定义元数据（如 FILTER、INFO、FORMAT 字段）
• 列头：#CHROM、POS、ID、REF、ALT、QUAL、FILTER、INFO、FORMAT、[样本名]
• 变异记录：每行描述一个变异位点

关键字段解析：

• REF/ALT：参考 / 变异碱基（ALT 可为多个，逗号分隔）
• QUAL：变异质量（Phred 值，Q=-10log₁₀(P)，P 为错误概率）
• INFO：变异注释（如 DP = 总深度，AF = 等位基因频率）
• FORMAT：样本数据格式（如 GT:AD:DP，分别表示基因型、等位基因深度、总深度）

2. 处理实操与代码模板

(1) Python：使用 PyVCF/cyvcf2 解析 VCF 文件

python

运行

# 安装：pip install PyVCF cyvcf2
import vcf

# 读取VCF文件
vcf_reader = vcf.Reader(open("variants.vcf", "r"))

# 统计SNP数量
snp_count = 0
for record in vcf_reader:
    if len(record.ALT[0]) == 1 and len(record.REF) == 1:  # 仅统计单碱基变异
        snp_count += 1
print(f"总SNP数: {snp_count}")

# 提取特定样本的杂合变异
sample_name = "Sample1"
het_variants = []
for record in vcf_reader:
    sample_data = record.genotype(sample_name)
    if sample_data["GT"] == "0/1":  # 杂合子
        het_variants.append(record)
print(f"{sample_name}的杂合变异数: {len(het_variants)}")

(2) 命令行工具：bcftools/vcftools 进行变异过滤与分析

bash

运行

# 安装：conda install bcftools vcftools

# 提取高质量变异（QUAL≥30且PASS过滤）
bcftools view -i "QUAL>30 & FILTER='PASS'" variants.vcf > high_qual.vcf

# 提取特定区域变异（chr2:10000-20000）
vcftools --vcf variants.vcf --bed regions.bed --out extracted_vars

# 计算样本等位基因频率
vcftools --vcf variants.vcf --freq --out allele_freq

四、GFF 格式：基因组注释的通用标准

1. 基本结构与原理

GFF（通用特征格式）是描述基因组特征（基因、外显子、启动子等）的标准格式，GFF3 版本由 9 列组成：

列	名称	描述
1	seqid	序列标识符（如染色体名）
2	source	特征来源（如数据库或预测软件）
3	type	特征类型（如 gene、exon、CDS）
4	start	起始位置（1-based）
5	end	终止位置
6	score	置信度得分（数字，可为 E 值 / P 值）
7	strand	链（+/-/.）
8	phase	仅用于 CDS，指示翻译相位
9	attributes	键值对属性（如 ID=gene001;Name=BRCA1）

2. 处理实操与代码模板

(1) Python：解析 GFF 文件提取基因信息

python

运行

# 安装：pip install biopython
from Bio import SeqIO
from Bio.SeqFeature import SeqFeature

# 读取GFF文件并提取所有基因特征
gene_features = []
with open("annotation.gff3", "r") as gff_file:
    for line in gff_file:
        if line.startswith("#"):  # 跳过注释行
            continue
        fields = line.strip().split("\t")
        if fields[2] == "gene":  # 仅提取基因特征
            # 解析属性
            attrs = {}
            for attr in fields[8].split(";"):
                key, value = attr.split("=")
                attrs[key] = value
            
            gene_features.append({
                "seqid": fields[0],
                "start": int(fields[3]),
                "end": int(fields[4]),
                "strand": fields[6],
                "id": attrs.get("ID"),
                "name": attrs.get("Name")
            })

print(f"共提取{len(gene_features)}个基因")

(2) 命令行工具：gffread 进行 GFF/GTF 格式转换与提取

bash

运行

# 安装：conda install gffread

# GFF转GTF
gffread annotation.gff -T -o output.gtf

# 提取特定类型特征（如CDS）
gffread annotation.gff -g genome.fa -x cds.fa -t CDS

# 统计各染色体基因数量
gffread annotation.gff -s | cut -f1,3 | sort | uniq -c

五、数据格式转换工作流

1. FASTQ → BAM：测序数据比对流程

plaintext

fastq → 质控(fastp) → 比对(bwa) → SAM → 排序/压缩(samtools) → BAM

核心命令：

bash

运行

# 质控
fastp -i input_R1.fq -I input_R2.fq -o clean_R1.fq -O clean_R2.fq

# 比对（以人类基因组为例）
bwa mem -t 8 hg38.fa clean_R1.fq clean_R2.fq > alignment.sam

# 转换为BAM并优化
samtools view -Sb alignment.sam | samtools sort -o sorted.bam
samtools index sorted.bam

2. BAM → VCF：变异检测流程

plaintext

BAM → 变异检测(samtools/bcftools) → 过滤 → VCF

核心命令：

bash

运行

# 变异检测
samtools mpileup -f hg38.fa sorted.bam | bcftools call -mv -o raw_variants.vcf

# 质量过滤
bcftools filter -i "QUAL>30 & DP>10 & AF>0.05" raw_variants.vcf -o filtered_vcf.vcf

六、实用代码模板汇总

1. FASTQ 质量控制（Python）

python

运行

# 基于Biopython的FASTQ质量过滤函数
def fastq_quality_filter(input_file, output_file, min_qual=20, min_len=30):
    """
    过滤FASTQ文件，保留质量值≥min_qual且长度≥min_len的reads
    """
    with open(input_file, "r") as in_handle, open(output_file, "w") as out_handle:
        for title, seq, qual in FastqGeneralIterator(in_handle):
            # 计算平均质量值
            avg_qual = sum(ord(c) - 33 for c in qual) / len(qual)
            if avg_qual >= min_qual and len(seq) >= min_len:
                out_handle.write(f"@{title}\n{seq}\n+\n{qual}\n")

# 使用示例
fastq_quality_filter("raw.fastq", "clean.fastq", min_qual=25, min_len=50)

2. BAM 覆盖率统计（Python）

python

运行

# 计算BAM文件在各染色体上的覆盖率
def bam_coverage_stats(bam_path):
    bamfile = pysam.AlignmentFile(bam_path, "rb")
    ref_names = bamfile.references
    
    coverage_stats = {}
    for ref_name in ref_names:
        # 计算该染色体的总覆盖碱基数
        total_covered = 0
        total_length = bamfile.lengths[bamfile.get_tid(ref_name)]
        
        # 计算每个位置的覆盖深度
        for pos, depth in enumerate(pysam.depth(bamfile, ref_name)):
            if depth > 0:
                total_covered += 1
        
        coverage_rate = total_covered / total_length * 100
        coverage_stats[ref_name] = {
            "total_length": total_length,
            "covered_bases": total_covered,
            "coverage_rate": coverage_rate
        }
    
    return coverage_stats

# 使用示例
stats = bam_coverage_stats("sorted.bam")
for chrom, stat in stats.items():
    print(f"{chrom}: 长度={stat['total_length']/1000:.1f}kb, 覆盖={stat['covered_bases']/1000:.1f}kb ({stat['coverage_rate']:.2f}%)")

3. VCF 变异注释（Python）

python

运行

# 从VCF文件中提取变异信息并添加功能注释
def vcf_annotation_parser(vcf_path, ref_genome_fa):
    """
    解析VCF文件，提取变异信息并预测功能影响
    """
    import numpy as np
    from Bio import SeqIO
    from Bio.Seq import Seq
    
    vcf_reader = vcf.Reader(open(vcf_path, "r"))
    ref_genome = SeqIO.read(ref_genome_fa, "fasta").seq
    
    variant_annotations = []
    
    for record in vcf_reader:
        chrom = record.CHROM
        pos = record.POS
        ref_base = record.REF
        alt_base = str(record.ALT[0])  # 仅考虑第一个变异
        
        # 提取变异前后的序列（以变异位点为中心，取10bp窗口）
        start = max(1, pos-10)
        end = pos+10
        ref_seq = ref_genome[start-1:end]  # 1-based to 0-based
        alt_seq = ref_seq[:pos-start] + alt_base + ref_seq[pos-start+1:]
        
        variant_annotations.append({
            "chrom": chrom,
            "pos": pos,
            "ref": ref_base,
            "alt": alt_base,
            "ref_seq": str(ref_seq),
            "alt_seq": str(alt_seq),
            "qual": record.QUAL,
            "filter": record.FILTER
        })
    
    return variant_annotations

# 使用示例
annotations = vcf_annotation_parser("variants.vcf", "hg38.fa")
for anno in annotations[:5]:  # 打印前5个变异
    print(f"变异位点: {anno['chrom']}:{anno['pos']} {anno['ref']}>{anno['alt']}")
    print(f"参考序列: {anno['ref_seq']}")
    print(f"变异序列: {anno['alt_seq']}\n")

七、总结与延伸

本手册涵盖了生信分析中四大核心数据格式的基本原理与实操方法：

• FASTQ：原始测序数据，包含序列和质量值，是所有分析的起点
• BAM：高效存储比对结果，支持快速检索和统计
• VCF：变异信息的标准格式，是基因组变异分析的核心
• GFF：基因组注释的通用标准，用于描述基因结构和功能