做全基因组测序WGS生信分析也有几年了从一开始拿到fastq文件手足无措到现在能完整跑通从质控到变异注释的全流程中间踩过的坑确实不少。这期间有不少朋友问我WGS的生信分析到底该怎么做从哪一步开始需要哪些工具参数怎么调。说实话网上的教程很多但要么是只讲某一步要么是流水账式地贴一遍代码真正能把原理、参数选择、踩坑经验串起来讲清楚的并不多。这篇就基于我自己的实操经验把WGS生信分析从原始数据到最终变异位点的完整链路拆开揉碎讲一遍希望能给刚入门或者正在被流程折磨的朋友一些参考。WGS也就是全基因组测序是基因组学研究里覆盖面最广、信息量最大的一种测序方式。它不像靶向测序那样只盯着某个区域而是把整个基因组的30亿个碱基对都测一遍。这意味着它适合做全基因组范围的变异检测包括SNP、InDel、CNV、SV等还能用于构建系统进化树、进行群体遗传分析、寻找疾病候选基因等。对于做生信的人来说WGS的数据量最大、计算资源要求最高、分析流程最完整把这套流程摸透了再回去做WES、RNA-seq都会轻松很多。我下面要讲的这套流程是建立在Illumina平台、双端150bp读长、30x测序深度的前提下这也是目前最主流、最通用的WGS数据形态。整个分析链路包括数据质控、比对、排序去重、碱基质量校正、变异检测、变异过滤、变异注释。每一步我都会讲清楚为什么要做、怎么做、参数怎么选以及实际跑数据时容易踩哪些坑。1. WGS分析先搞清楚这几点原理与数据形态1.1 全基因组测序到底在测什么全基因组测序的核心逻辑是把基因组DNA打成小片段加上接头后上机测序产出大量短读长读长也就是read一般150bp然后通过生信手段把这些读长比对回参考基因组上找到每个位置上的碱基和参考基因组是否一致不一致的地方就是变异位点。这里有个关键概念要理解测序不是直接把整个基因组从头到尾读一遍而是把基因组随机打断成小片段然后对这些小片段的两端pair-end进行测序。打断是随机的所以一个基因组区域会被很多不同的片段覆盖每条read从不同的起始位置开始读这样比对回去的时候同一个位置有多条read覆盖才能保证碱基判定的准确性。这就引出了测序深度的概念。深度depth或coverage指的是基因组上每个位置被测到的平均次数。30x就表示平均每个位点被30条read覆盖。这个数字不是随便定的它背后的逻辑是碱基判定是一个概率问题覆盖度越高错误判定的概率越低。在30x深度下一个纯合变异的位点理论上会有接近100%的read显示变异杂合变异大约有50%的read显示变异这些比例关系是后续变异检测可信度的基础。我刚才说的打断、测序、比对、找变异这个逻辑链是理解WGS生信分析的骨架。后面每一步操作都是在为这个逻辑链上的某个环节服务。质控是为了保证从测序仪出来的数据干净可靠比对是为了把read定位回基因组坐标系变异检测是把可靠的差异位点找出来注释是搞清楚这些差异在不同基因组元件中的位置和可能的生物学影响。1.2 一个WGS样本的数据量有多大在跑通流程之前你最好对数据量有个心理预期。一个人全基因组的长度大约3.1Gb30x深度意味着测序仪要产出大约90-100Gb的原始碱基数据。换算成fastq文件一个样本的双端reads在4到5亿条左右fastq文件大约占200GB的存储空间。我最早跑WGS的时候对数据量完全没概念以为一两个小时就能跑完结果光比对这一步就跑了一个通宵。后来才总结出经验一个30x的WGS样本从fastq开始到最终得到注释好的vcf文件在32核服务器上大约需要40到60个小时如果机器配置差一些或者磁盘IO不够快耗时会更长。这个时间估算很重要能帮你在安排项目进度的时候做合理规划。计算方面比对和变异检测是两个最吃资源、最耗时的步骤。比对建议至少分配16到32个线程内存不低于32GB变异检测阶段GATK HaplotypeCaller对内存的需求更高建议64GB以上否则很容易OOM内存溢出这一点我后面还会提到。2. 生信分析流程设计与工具选型2.1 标准化流程的骨架从fastq到vcf现在跑WGS分析基本上都会走GATK官方推荐的Best Practices流程这一套流程被验证过无数次可靠性有保障。整套流程可以拆成三段数据准备段、变异发现段、变异整理段。数据准备段的输入是原始测序fastq文件测序仪直接下机的数据输出是经过清洗、排序、去重、质量校正的比对文件BAM文件。这个阶段的核心目的是把reads准确可靠地放到基因组上并消除技术误差。变异发现段的输入是BAM文件输出是包含每个位点基因型信息的gVCF文件。这个阶段的核心是HaplotypeCaller它会根据有向图算法重新组装每个区域的单倍型再进行碱基判定这样能显著提高InDel检测的准确性。变异整理段的输入是所有样本的gVCF文件输出是最终的VCF文件。这个阶段包括联合基因分型GenomicsDBImport GenotypeGVCFs和VQSRVariant Quality Score Recalibration质量控制目的是把多个样本的数据合并起来做统一的基因型判定和质量校正。这条流程的好处是可以直接从官方拿到最新更新和完整的参数文档社区用户多遇到问题基本都能找到解决方案。我个人建议新手就用这套标准流程不要一上来就去折腾那些所谓的“简化版”“加速版”流程等你把标准流程跑通了了解每一步的输入输出是什么、参数为什么这么设再去优化也不迟。2.2 核心工具选型为什么是BWA GATK比对工具目前主流的选择是BWA-MEM也包括它对较长读长的变体BWA-MEM2。BWA-MEM2利用SIMD向量化技术加速了种子搜索和扩展过程在保持与BWA-MEM几乎一致比对精度的前提下速度接近两倍我后来一直用的是BWA-MEM2。为什么选BWA而不是其他工具因为它对Illumina短读长的比对速度和准确率是均衡得最好的。而且在WGS场景中我们不仅要比对还要考虑后续GATK的兼容性BWA输出的SAM/BAM文件可以被GATK直接读取不需要额外转换这一条就省了很多事。变异检测这边GATK是绕不开的。虽然也有freebayes、samtools mpileupbcftools等替代方案但GATK的优势在于它有一套完整且统一的质量控制体系。VQSR可以用已知位点数据库对变异位点的质量打分进行校正这在大型WGS项目中能显著降低假阳性率。如果你的项目样本量够大建议至少30个WGS样本VQSR的收益很明显样本数不够时用硬过滤Hard Filter也够用GATK官方已经给出了推荐的过滤参数。我这里补充一个工具选型的逻辑在做生信流程设计的时候不要只盯着某个工具单步跑得有多快更要看工具之间的衔接是否顺畅、参数体系是否统一、出了问题好不好排查。BWA GATK这套组合经过大量项目验证各个步骤之间的接口都很干净出了问题能很快定位到具体环节。2.3 参考基因组的选择与准备参考基因组是比对的标准坐标系选择不当会直接影响变异检测结果。目前人类基因组最常用的是GRCh38也就是hg38相比之前的GRCh37hg19主要更新了几千个错配位点和一些结构变异的表示方式。如果是从零开始的WGS项目建议直接用GRCh38GATK的已知位点数据库、人群频率数据库gnomAD等现在都默认基于GRCh38。参考基因组准备这件事看起来简单但其实有几个细节需要注意。第一个细节是一定要去GATK官方的资源包页面下载b37/hg38对应的参考基因组文件不要随便从网上某个地方下载一个不知道什么版本的fasta。参考基因组的版本不对后面所有比对结果都是错的而且这种错误非常隐蔽不容易发现。第二个细节是参考基因组必须建立索引文件包括.fai索引和.dict字典文件如果是BWA比对还要提前建好BWA索引。这些索引文件是一次性建立的但如果没有它们比对和变异检测工具会直接报错。另外建议给参考基因组建一个专门的目录把fasta、索引文件、已知变异数据库都放一起。我见过不少人把参考基因组文件散落在各个项目目录里每次跑新项目都要重新指定路径既容易出错又不方便统一管理。我这里习惯建一个reference目录里面有fasta、dict、fai、bwa索引以及dbsnp、Mills、1000G等已知位点的vcf.gz文件以后的每个项目都用这个目录省心很多。3. 实操过程与关键环节实现3.1 质控什么样的数据能进入分析流程质控是WGS分析的第一步也是最不能省的一步。测序仪下机的数据质量再高也会存在低质量碱基、接头序列污染、重复序列等问题如果不过滤掉后面比对和变异检测的结果都会被影响。常用的质控工具是FastQC和MultiQC。FastQC会对输入的fastq文件做一套全面的质量评估包括每个位置的平均碱基质量、GC含量分布、接头污染情况、重复率、Kmer富集等。MultiQC则能把多个样本的FastQC报告汇总成一个交互式HTML方便对比和筛选。拿到质控报告之后我一般会重点看几个指标。一是每个位置的碱基质量分布NGS测序质量会随读长增加而下降正常情况前100bp的Q30比例应该在90%以上后面50bp略低是正常的。二是GC含量分布人类基因组的GC含量大约是41%如果实测GC含量出现多峰分布提示可能有不同物种的DNA污染。三是接头含量现在很多测序项目在交付前已经做了去接头处理但如果原始数据里接头含量超过1%就需要用cutadapt或Trimmomatic处理。过滤这一步我的建议是宁缺毋滥。对于双端reads如果某一条read的质量太差可以直接把这一对read都丢弃因为比对算法需要双端信息来确定插入片段大小只保留单端数据反而会引入噪音。另外不要一刀切地把所有Q值低于某个阈值的碱基都切掉这样做会过度缩短reads反而降低比对率。我一般只做两端柔软切除soft clipping和整条read的质量筛选中段质量略低不影响比对结果。如果你的测序数据是别人交付的记得问清楚测序平台、读长、插入片段大小、测序深度这几个关键参数。拿到数据后先做一轮FastQC判断数据形态是否符合预期再决定是否直接进入下一步。我之前接过一批从测序公司拿到的WGS数据对方声称是PE150结果FastQC跑完发现reads长度只有128bp且3端质量持续降低原来是测序过程中出现了异常终止这批数据最后只能降级使用损失不小。3.2 比对环节BWA-MEM2参数与执行细节比对这一步做的是将质控后的reads比对回参考基因组得到SAM文件然后转换为BAM文件。BWA比对的核心命令大概是这样的形式# 建立BWA索引如果还没建的话 bwa-mem2 index GRCh38.fa # 执行比对-t指定线程数 bwa-mem2 mem -t 32 -R RG\tID:sample1\tSM:sample1\tPL:ILLUMINA\tLB:sample1 GRCh38.fa sample1_R1.fastq.gz sample1_R2.fastq.gz sample1.sam这里有个极其重要的细节-R参数一定要写完整。Read Group信息里包含了样本ID、文库ID、平台类型等信息GATK在后续的重复标记、碱基质量校正、变异检测中都会用到这些标签。如果你不写Read GroupGATK会默认给你一个很随意的值会导致后续步骤报错或者样本混淆这个错我踩过一次后来凡是写脚本都会把-R参数单独列出来检查无误再提交任务。比对完成后得到的SAM文件需要转换成BAM格式SAM是文本格式体积庞大BAM是二进制压缩格式处理速度更快然后还要做排序因为比对输出的reads顺序是乱的而GATK要求reads按坐标排序。这一步推荐用samtools完成# 转换并排序 samtools view -bS - 16 sample1.sam sample1.bam samtools sort - 16 -m 4G -o sample1.sorted.bam sample1.bam samtools index sample1.sorted.bam排序的参数可以关注一下。-m参数控制每个线程最多用多少内存32个线程的时候就给了4G这样总内存上限是128G在常见的高配置服务器上没问题。如果机器内存紧张可以把-m调低一些但排序速度会相应变慢。比对完成后用一个叫samtools flagstat的常用命令来检查比对结果的质量看比对率是否正常。对于人类WGS数据比对率mapped ratio通常在99%以上。如果比对率明显偏低低于95%说明参考基因组版本可能不对、样本可能不是人类的、或者reads质量太差需要回头排查不可大意。3.3 去重与碱基质量校正比对这一步做完BAM文件里还有一类技术噪音需要处理就是PCR重复PCR duplicates。在文库构建过程中DNA片段经过PCR扩增同一个原始片段可能会被扩增出很多份测序后这些reads会被比对回同一个位置。如果不去除重复reads变异检测时这些重复reads的大量投票会让某些位点的覆盖度虚高影响等位基因频率计算的准确性。去除重复reads用GATK的MarkDuplicates工具它通过比对位置和UMI唯一分子标签来识别重复。对人类WGS项目我处理过的数据里重复率通常在5%到20%之间。如果重复率超过30%需要留意文库质量或扩增条件是否出了问题。这一步完成后可以顺手用samtools flagstat再确认一下有效reads的比例。紧接着是BQSRBase Quality Score Recalibration碱基质量分数重新校正。这个步骤是GATK流程里比较有特色的一个环节目的是利用已知变异位点数据库消除测序仪产生的系统性碱基质量误差。测序仪的碱基质量分数是基于测序仪内部算法计算出来的会受批次效应、试剂、跑胶条件等因素影响导致质量分数存在偏差。BQSR会重新校准这些分数让它们更接近真实错误概率。BQSR的过程分为两步先是用BaseRecalibrator根据已知位点数据库和协变量信息计算出校正表再用ApplyBQSR把校正表应用到BAM文件上。这步操作对变异检测的准确性提升在WGS数据上的收益比WES更明显因为WGS覆盖的区域更广校正能应用到更多位点上。BQSR的计算速度也较快但会输出一个新的BAM文件注意给磁盘留足空间。3.4 变异检测HaplotypeCaller的正确打开方式变异检测是WGS分析的核心环节也是整个流程中最耗时的步骤。GATK的HaplotypeCaller采用的方法是对每个区域先用de Bruijn-like图模型进行局部组装生成可能的单倍型再把每条read比对回这些单倍型上最后基于这些单倍型的支持情况计算每个位点的基因型。HaplotypeCaller有两种推荐运行模式。一种是单样本直接生成VCF的模式适合快速看结果另一种是生成gVCF的模式先对每个样本单独运行再用GenomicsDBImport合并最后用GenotypeGVCFs联合基因分型。推荐后者因为生成gVCF之后如果后续有新的样本加入只需要对新样本跑HaplotypeCaller再合并一次就行不需要回到最初重新计算所有样本。HaplotypeCaller是公认的计算密集型任务建议分配足够的资源gatk --java-options -Xmx64G HaplotypeCaller \ -R GRCh38.fa \ -I sample1.sorted.dedup.bqsr.bam \ -O sample1.g.vcf.gz \ -ERC GVCF \ --native-pair-hmm-threads 16这里有两个参数值得展开讲一下。一个是--native-pair-hmm-threads它控制PairHMM模型的线程数这个环节是HaplotypeCaller最耗时的部分通过将它设为较高的值比如16能显著缩短运行时间。另一个是-Java-Xmx参数HaplotypeCaller对Java堆内存有一定需求设小了会直接OOM但设太大会导致内存浪费一般设成可用物理内存的一半左右即可。如果你有多台机器可以用建议在不同机器上并行跑不同染色体的gVCF最后再合并这样能进一步缩短整体时间。GATK官方支持按区间interval切分任务最常见的做法是按染色体号分别运行1号染色体单独跑22号染色体等可以直接合并几条一起跑以实现计算资源的均衡利用。3.5 联合基因分型与硬过滤当所有样本的gVCF都生成之后接下来就是联合基因分型。先导入GenomicsDB# 建立样本映射文件格式为 样本名 gvcf路径 gatk --java-options -Xmx64G GenomicsDBImport \ -V sample_map.tsv \ --genomicsdb-workspace-path my_database \ -L intervals.list \ --reader-threads 5这一步对内存要求很高需要根据数据库规模调整Java堆内存我跑过32个样本的联合导入64G内存基本够用。如果样本量很大建议按染色体分批导入避免单次任务内存溢出。导入完成后gatk --java-options -Xmx32G GenotypeGVCFs \ -R GRCh38.fa \ -V gendb://my_database \ -O cohort.g.vcf.gz联合基因分型完成后得到的VCF文件还需要做质量控制过滤掉那些可能是假阳性的变异位点。质量控制有两种方式VQSR和硬过滤。VQSR适合样本量大的项目官方建议至少30个WGS样本利用已知位点数据库训练一个高斯混合模型给每个位点的质量分数打分样本量不足时硬过滤更可靠。硬过滤常用参数可以参考GATK官网推荐的配置文件核心是QD、FS、MQ、MQRankSum、ReadPosRankSum这几个指标。比如SNP的推荐阈值是QD2.0、FS60.0、MQ40.0、MQRankSum-12.5、ReadPosRankSum-8.0。这些指标分别反映了变异的测序质量、链偏倚、比对质量、映射质量与邻位点的差异等每个指标都有生物学意义不是随便设的数字。如果你用硬过滤建议对照GATK官方最新推荐的参数不要用过时的数值。3.6 变异注释从位点到功能变异检测输出的VCF文件里是一堆位点和基因型信息要理解这些位点的生物学意义还需要对它们进行功能注释。功能注释要做的事包括判断变异位点是否落在基因区域内是外显子还是内含子是否是错义突变、无义突变、剪接位点变异等以及这个位点在人群中的频率是多少。目前最常用的注释工具是ANNOVAR和VEPVariant Effect Predictor。ANNOVAR的命令简洁高效例如对GRCh38的人全基因组注释的核心命令是table_annovar.pl cohort.vcf humandb/ \ -buildver hg38 \ -out cohort_anno \ -remove \ -protocol refGene,1000g2015aug_all,clinvar_20221231,gnomad312_genome \ -operation g,f,f,f \ -nastring . \ -vcfinput这里的-protocol参数列出了参考基因集和各个数据库的名称-operation参数标明每个数据库的类型g表示基因相关f表示频率相关。我自己一般会把refGene、gnomAD、ClinVar、1000G这几个都加上频率数据库能帮忙筛掉那些人群高频的普通多态位点ClinVar能直接标注出已知的致病或良性位点。除了这些基础注释如果你是做疾病研究可能还需要额外注释OMIM、GWAS Catalog、dbSNP等数据库如果你是做肿瘤WGS可以加上肿瘤相关的数据库如COSMIC、ICGC等以及可以用MSI、TMB等指标来评估肿瘤突变负荷。注释这一步的输出是一个带有多列注释信息的VCF或者表格文件之后的统计分析和筛选都基于这个文件展开。4. 常见问题与排查技巧实录4.1 比对率低或者GC偏差异常如果比对后比对率明显低于预期或者FastQC报告显示GC含量分布不合理需要先排查参考基因组版本是否正确、数据是否有交叉污染。交叉污染在WGS数据中的表现是杂合位点数异常偏高可以用VerifyBamID来检查。我自己遇到过一次比较典型的情况某个样本的比对率只有96%反复排查后发现是参考基因组版本和测序公司使用的版本不一致公司用GRCh37我用了GRCh38重新比对了正确的版本之后比对率回到了99.3%。所以拿到数据之前一定要先确认参考基因组版本。4.2 HaplotypeCaller内存溢出或运行卡死HaplotypeCaller是流程中最容易出问题的环节主要原因有两个内存不足和某个区域read堆积过载。内存不足直接报OOM这时需要调大-Xmx参数或者把任务按染色体切分后再跑。区域read堆积过载通常表现为某个区间运行速度极慢比如chr1上某个高重复区域。这类区域本身比对到的reads就非常多HaplotypeCaller在做局部组装时会消耗大量内存和时间。解决办法是给HaplotypeCaller传入一个排除区间列表比如已知的高重复区域跳过这些区域或者把该区域单独拆出来跑用更大的内存应对。还有一个容易忽略的点HaplotypeCaller对输入BAM文件有要求sample的性别、是否混样等信息要标注清楚。如果你用-gender参数指定了性别但BAM头部的SM标签不匹配HaplotypeCaller可能不会报错但结果会受影响尤其是X和Y染色体的变异检测。4.3 硬过滤后变异位点数量异常过滤后SNP数量明显偏少或偏多说明过滤阈值或者数据质量有问题。正常情况下一个人全基因组的SNP数量在350万到450万左右InDel在50万左右这个范围可以作为参考。如果SNP数量远多于400万要怀疑是否有DNA污染两个样本混在一起了或参考基因组选择错误如果远少于350万可能是BAM文件深度不足或者该样本与参考基因组亲缘关系太近。我跑过一个祖源研究项目样本来自同一家族他们之间的SNP数量会比正常不相关个体少一些这是正常现象不属于异常。4.4 不同批次WGS数据的合并分析问题WGS项目经常是分批做测序的不同批次的文库构建、测序条件可能不同这种批次效应会在下游分析中引入假阳性信号。处理方式是在数据合并前做质量控制尽量保证各批次的质量指标相近在分析阶段可以用主成分分析或聚类分析来检查样本是否存在明显的批次分层如果要做关联分析可以将批次作为协变量纳入模型。另外不同批次的数据在变异检测阶段要统一使用同一个参考基因组版本、同一套过滤参数、同一个数据库版本。我见过有人在两批数据中用了不同版本的ClinVar数据库结果注释出来的致病基因列表对不上排查了很久才发现是数据库版本不一致。4.5 数据存储和备份的管理经验最后说一下存储问题这个在WGS项目中经常被忽略。一个30x WGS样本的fastq文件约200GB比对排序后的BAM约100GB去重和BQSR后的BAM又是100GBgVCF大约1-2GB最终注释好的VCF几十到几百MB。一个样本算下来原始中间文件的总存储需求能到400-500GB20个样本就是接近10TB。我的建议是中间文件可以按需保留但fastq和最终的VCF一定要备份。fastq是原始数据无法重新生成丢了就只能重新测序VCF是分析成果后续做任何统计和解读都靠它。而中间产生的BAM文件如果计算资源和时间允许建议至少保留去重BQSR后的BAM因为它版本可重复利用后续换一套GATK版本或加一个新的变异检测工具直接从BAM开始不用重新比对。存储方面可以考虑使用压缩文件系统、重复数据删除技术或者把不常用的中间文件迁移到冷存储上这样能省不少成本。另外每个样本的Read Group信息和项目参数我都建议记在一个元数据表格里这看起来是小事但一旦后面要继续分析或者跟别人协作这份元数据就是救命稻草。5. 写在最后的一点个人经验全基因组测序的生信分析本质上是一个数据加工流水线。从测序仪出来的原始数据经过质控、比对、变异检测、注释最终变成一张张含有生物信息的表格。这条流水线的每一环都需要严谨对待因为上游一个小小的问题到下游可能就被放大成完全错误的研究结论。我个人的体会是跑通流程只是第一步真正有价值的收获来自于对每一步原理的理解和对异常结果的敏感度。多花点时间看看GATK官方文档多把自己的数据和别人的结果对比多问几个为什么这套流程就会越用越顺手。最后再分享一个小技巧不管项目多急我在正式跑全量数据前一定会先用一条染色体的小数据集把整个流程完整走一遍确认每一步的输入输出都没问题、资源占用合理再提交全量任务。这样做一次全流程测试能省掉无数个排错和重启任务的深夜。