全基因组重测序上游分析流程--随笔15-洪萨配资

全基因组重测序上游分析流程｜从软件部署到变异检测，超细致实操指南

作为科研新手，第一次上手全基因组重测序数据处理时，我踩过不少软件安装的坑、碰过参数设置的雷。如今整理出这份超详细流程，从前期准备到最终变异过滤，每一步都标注了关键注意事项，跟着练一遍就能快速上手。觉得有用的话，别忘了点赞收藏哦～

适用场景：动植物全基因组重测序上游分析（变异检测核心流程）

核心工具：BWA、Samtools、Picard、GATK4（全开源，附conda安装命令）

0 前期准备：软件部署与环境搭建

重测序上游分析的核心工具就4个，全部开源可通过conda快速安装，推荐Linux集群或Mac OS系统，Windows建议用WSL2。重点注意GATK4的版本特性——虽然是新版本，但100%开源且适配大规模数据，是未来主流，本文全程基于GATK4构建流程。

工具名称	核心功能	conda安装命令	注意事项
BWA	NGS数据与参考基因组比对	conda install -c bioconda bwa	C语言编写，需系统支持编译
Samtools	BAM/SAM文件处理（排序、索引等）	conda install -c bioconda samtools	与Picard功能互补，必装工具
Picard	标记重复序列、文件格式处理	conda install -c bioconda picard	Java编写，需Java 1.8+环境
GATK4	变异检测、基因型推断	conda install -c bioconda gatk4	4.x比3.x更适配集群，全开源

避坑指南：安装生信软件务必加-c bioconda指定频道，避免下载到旧版本或错误包。Java版本过低会导致Picard报错，可通过java -version检查，低于1.8则需重新安装。

1 原始数据质控：快速过一遍，心里有底

现在测序公司交付的基本都是经过初步处理的clean data，但自己验证一步更放心。核心用两个工具：

FastQC：可视化展示数据质量（碱基分布、测序错误率等），命令：fastqc read_1.fq.gz read_2.fq.gz
fastp：批量清洗数据（去除接头、低质量碱基），命令：fastp -i read_1.fq.gz -I read_2.fq.gz -o clean_1.fq.gz -O clean_2.fq.gz

如果FastQC报告中“Per base sequence quality”出现红色区域，或接头污染率高，就需要用fastp调整参数（如--cut_front去除前端低质量碱基）再处理。

2 核心流程：从序列比对到变异检测

这部分是重测序分析的核心，每一步都有明确的逻辑和避坑点，跟着步骤走准没错。

2.1 序列比对：给短读长“找家”

NGS测出来的短序列（read）是随机打乱的，必须通过比对找到它们在参考基因组上的位置。BWA是目前最权威的工具，核心靠“索引构建+比对”两步。

步骤1：构建参考基因组索引

索引能让BWA快速定位序列，相当于给参考基因组建“目录”。命令超简单：

bwa index genome.fasta

运行后会生成5个以genome.fasta为前缀的文件（.amb/.ann/.bwt等），这些是比对的关键，别删！

步骤2：双末端序列比对

重测序常用双末端测序（PE），两个fq文件分别对应DNA片段的两端，比对时要一起输入。这里有个超级关键的参数——Read Group（-R），直接影响后续GATK分析！

bwa mem -t 4 -R '@RG\tID:lane1\tPL:illumina\tLB:lib1\tSM:sample1' genome.fasta clean_1.fq.gz clean_2.fq.gz | samtools view -S -b - > sample1.bam

ID：测序lane编号（从fq文件名获取，如lane1）
PL：测序平台（必须是GATK认可的，如illumina、COMPLETE，不能写“CG”“MGI”）
SM：样本ID（唯一标识，多样本分析时必用）
LB：文库名（可选，从测序报告获取）

避坑重点：-R参数的4个核心信息平台写错会报“not a recognized platform”错误，后期改起来很麻烦！

命令解析：-t 4用4个线程加速；管道符|直接将比对结果（SAM格式）转给Samtools，用-b转为二进制BAM格式（节省空间，后续分析更高效）。

2.2 数据排序：按染色体位置“排好队”

BWA比对后的BAM文件是按read的测序顺序排列的，而后续分析需要按染色体位置排序。用Samtools完成，命令：

samtools sort -@ 4 -m 4G -O bam -o sample1.sorted.bam sample1.bam

参数说明：-m 4G限制每个线程用4G内存，避免服务器内存溢出；文件名加“sorted”标识，后续好区分。排序后文件会略小，是压缩算法导致的，内容无损失。

2.3 标记重复序列：剔除PCR扩增的“赝品”

建库时的PCR扩增会产生大量重复序列，这些序列会干扰变异检测（增大假阳/假阴率），必须标记或去除。主流用Picard的MarkDuplicates，默认只标记不删除（更灵活）。

picard MarkDuplicates I=sample1.sorted.bam O=sample1.sorted.markdup.bam M=sample1.markdup_metrics.txt

参数说明：I是输入文件，O是输出文件，M是重复序列统计报告（可查看重复率，一般低于30%算正常）。如果非要删除重复序列，加REMOVE_DUPLICATES=true参数即可。

2.4 构建索引：让工具“随机访问”文件

标记重复后的BAM文件需要建索引，方便后续工具快速定位特定区域。同时要给参考基因组做GATK专用索引，两步命令：

# 给BAM文件建索引 samtools index sample1.sorted.markdup.bam # 给参考基因组建GATK索引（生成.dict和.fai文件） gatk CreateSequenceDictionary -R genome.fasta -O genome.dict && samtools faidx genome.fasta

运行后会生成sample1.sorted.markdup.bam.bai（BAM索引）、genome.dict和genome.fasta.fai（参考基因组索引），这三个文件缺一不可。

2.5 变异检测：从GVCF到最终VCF

GATK的HaplotypeCaller是目前最优的变异检测工具，支持单样本和多样本分析，核心分“生成GVCF→合并→基因型推断”三步。

步骤1：单样本生成GVCF

GVCF文件包含所有位点信息（无论是否变异），便于后续多样本合并分析。命令：

gatk HaplotypeCaller -R genome.fasta -I sample1.sorted.markdup.bam --emit-ref-confidence GVCF --min-base-quality-score 10 -O sample1.chr1.g.vcf.gz

如果样本多、染色体多，建议写shell脚本批量运行（循环修改染色体号和样本名），效率翻倍。

步骤2：合并多样本GVCF

多个样本按染色体合并，先把同染色体的GVCF文件名存成列表，再用CombineGVCFs合并：

# 生成GVCF列表 ls *.chr1.g.vcf.gz > chr1_gvcf.list # 合并 gatk CombineGVCFs -R genome.fasta -V chr1_gvcf.list -L 1 -O chr1.merged.g.vcf.gz

步骤3：基因型推断（生成VCF）

将合并后的GVCF转为最终的变异文件（VCF），包含SNP和InDel信息：

gatk GenotypeGVCFs -R genome.fasta -V chr1.merged.g.vcf.gz -O chr1.genotype.vcf.gz

2.6 变异过滤：剔除假阳性，保留可靠结果

刚生成的VCF是“原始数据”，包含大量假阳性变异，需要过滤。分SNP和InDel两类处理，非人类物种建议用“硬过滤”（人类可用VQSR，依赖已知变异集）。

# 提取SNP gatk SelectVariants -R genome.fasta -V chr1.genotype.vcf.gz -O chr1.snp.vcf -select-type SNP # 过滤SNP（核心参数，可根据数据调整） gatk VariantFiltration -V chr1.snp.vcf -O chr1.snp.filter.vcf -R genome.fasta \ --filter-expression "QD < 2.0 || FS > 60.0 || MQ < 40.0" \ --filter-name "SNP_filter"

过滤参数说明：QD（变异质量值）、FS（碱基偏倚）、MQ（比对质量），这些是GATK推荐的核心指标，过滤后标记为“SNP_filter”的位点就是需要剔除的假阳性。