1. 从"有没有变异"到"变异在哪里":WGS的核心定位
做了几年生信,被问得最多的一个问题就是:WGS到底比靶向测序(比如全外显子组测序WES、Amplicon panel)强在哪?很多刚接触测序数据的同学,把WGS等同于"测得更深、更贵、数据量更大",这个理解不能说错,但没有触及WGS真正改变分析思路的地方——WGS是唯一一种以全基因组为对象、不依赖探针捕获、能够同时覆盖编码区和非编码区的测序策略。这句话听起来像概念搬运,实际落地之后你会发现它改变了整个分析流程的设计逻辑。
先说说WGS能干什么。单核苷酸变异(SNV)和小片段插入缺失(InDel)是最常见的诉求,这跟WES重叠;但WGS还有三个典型的差异化价值:结构变异(SV)的检测灵敏度远超WES,因为SV的断裂点往往落在内含子或基因间区,探针捕获根本拉不到;拷贝数变异(CNV)的断点分辨率可以做到更精细,尤其对肿瘤样本,WGS数据能直接推导出等位基因特异性拷贝数(allele-specific copy number),这是WES很难实现的;另外就是非编码区的调控变异,这类位点以前被当作"垃圾",现在越来越多的功能研究和临床解读开始关注它们。
所以,如果你手里已经有WGS数据,接下来这篇内容会完整走一遍从原始数据到变异解读的分析链路;如果你还在纠结要不要上WGS,我也会给出选型层面的判断依据。这里默认覆盖度为30x左右的标准人全基因组测序数据,这是目前科研和临床最常用的配置。
2. 建库、上机到数据产出:真正影响结果质量的环节
很多教程直接从fastq开始讲,似乎测序数据是天上掉下来的。但根据我这些年经手的项目来看,至少30%的分析问题和上游建库上机直接相关——接头污染、duplication偏高、碱基质量掉底、GC偏好性,这些在QC阶段才暴露的问题,根子都在实验环节。生信人员至少要理解上游流程,才知道哪些问题能通过分析手段修正,哪些必须返回去重新测序。
2.1 Illumina与华大平台的建库差异
主流的WGS建库思路都是:DNA片段化 → 末端修复 → 加A尾 → 连接接头 → PCR扩增(可选) → 上机测序。但Illumina和华大DNBSEQ两个平台在关键技术路线上有本质区别。
Illumina用的是桥式扩增(bridge amplification),DNA片段两端连接Y型接头后,通过flow cell表面固定的寡核苷酸链进行局部扩增,形成DNA簇(cluster)。DNBSEQ则是把连接好接头的环化DNA通过滚环扩增(RCA)生成DNA纳米球(DNB),吸附到阵列芯片上。DBN技术有一个天然优势:每个DNB只占一个位点,信号点是单拷贝,测序过程不会累积错误;Illumina的cluster是同步扩增的,虽然信号强度更高,但扩增错误会堆叠。
这个差异直接反映在数据上:DNBSEQ平台的duplication rate普遍更低,且对高GC区域有更好的覆盖均匀度。工业界的经验值是,同一样本用DNBSEQ跑出来的有效数据比例,通常比Illumina高出几个百分点。
2.2 PCR-free建库带来的数据变化
PCR-free(无PCR扩增)是当前WGS建库的主流做法。它跳过PCR步骤直接上机,保留原始DNA分子的多样性。传统建库中的PCR有几个副作用:高GC区域扩增效率低导致覆盖偏倚,重复序列被指数级放大的概率变高,PCR聚合酶的碱基错误被复制进最终数据。做全基因组测序,理论上不需要扩增这一步来达到足够的信号强度,所以PCR-free在WGS项目里几乎成了默认选项。
需要注意的是,PCR-free对起始DNA质量和总量要求更高。一般建议总量500ng以上、OD260/280在1.8-2.0之间、无明显的RNA污染。如果你的样本是FFPE或者微量DNA,做不了PCR-free,就只能走有PCR的流程。这时候,后续分析中的MarkDuplicate步骤就变得格外重要。
2.3 测序深度、读长与数据量的换算
30x WGS的标准含义是,人类基因组约3.1Gb,30x意味着需要约93Gb的有效数据(mappable bases)。但注意"有效"两个字——实际获得的总碱基数会高于这个值,因为存在read比对不上、duplicate、quality cutoff之后被过滤掉的部分。
平台选型层面有一个比较实用的对照表:
| 平台 | 读长 | 单端reads量(约) | 30x所需数据量 | 典型运行时间 |
|---|---|---|---|---|
| Illumina NovaSeq X | PE150 | 可灵活配置 | 约1Tb | 约2天 |
| DNBSEQ-T7 | PE150 | 可灵活配置 | 约1Tb | 约1-2天 |
| MGI DNBSEQ-G400 | PE150 | 中等通量 | 约300-600Gb | 约2天 |
| PacBio HiFi | 15-25kb长读 | 覆盖度更密集 | 约20x即可 | 数小时-1天 |
如果你有30x左右的WGS fastq数据,通常的文件体积是:fastq.gz约60-80GB(每个样本)。PE150是当前WGS最标准的读长配置,PE100也能做,但后续处理SV和短片段InDel时能明显感觉到差异。
3. 从fastq到BAM的"三级跳":比对、排序去重、碱基质量校正
拿到下机数据之后,真正的生信分析流程从这里开始。整个数据预处理的核心目标是把测序read"对齐"到参考基因组上,标记出质量可疑的重复片段,并校正系统性的碱基质量误差。每一步的产出都会作为下一步的输入,所以这部分的稳妥程度决定了后面变异检测的可靠性。
3.1 第一步:原始数据QC
用fastp或FastQC+MultiQC组合对fastq做质控,重点看几个指标:每条read的平均碱基质量(Q30百分比)、GC含量分布、接头残留率、重复率、插入片段长度分布。对于WGS数据,Q30比例一般应达到85%以上。如果接头残留超过5%,建议用fastp的先截后滤策略做一次轻量清洗。
值得提醒的是,现在的WGS测序服务商交付的数据通常已经做过demultiplex和adapter trimming,但依然是"通常",不能默认。我自己有个习惯:不管供应商说处理过没有,拿到数据永远先跑一遍fastp确认接头情况。
3.2 第二步:BWA-MEM比对
比对是决定后续变异检出精度的关键一步。WGS短读长比对的标准工具是BWA-MEM(或它的更新版BWA-MEM2),它将每个read与参考基因组比对,输出SAM/BAM格式的比对结果。
命令大致是:
# 建立索引(人基因组约3.1Gb,需要约5-10分钟) bwa index -a bwtsw GRCh38.fa # 比对并转为BAM bwa mem -t 32 -R '@RG\tID:sample1\tSM:sample1\tLB:lib1\tPL:illumina' GRCh38.fa sample1_R1.fastq.gz sample1_R2.fastq.gz \ | samtools view -bS - > sample1.raw.bam-R参数指定read group,这一段很多人会忽略,但后续GATK流程判断样本身份和合并多lane数据时都依赖它,绝对不能省。如果同一份DNA在多个lane/多个flow cell上测了,这一步不写清楚,后面合并时会有大麻烦。
人全基因组的比对速度大约为单线程每小时几百Mb到1Gb,实际取决于机器CPU。用32线程跑一份30x WGS,通常需要40分钟到1.5小时。比对阶段是最耗CPU的步骤之一,建议用高性能计算节点或者云上CPU密集实例。
3.3 第三步:排序与MarkDuplicate
比对后的BAM文件需要按参考基因组坐标排序,这一步用samtools sort。然后标记PCR或光学重复(duplicate reads)。WGS项目中duplication率一般在5%-15%之间,具体受建库质量和测序Cluster密度影响。
samtools sort -@ 16 -m 2G sample1.raw.bam -o sample1.sorted.bam # 标记重复 picard MarkDuplicates \ I=sample1.sorted.bam \ O=sample1.markdup.bam \ M=sample1.markdup_metrics.txt \ REMOVE_DUPLICATES=false关于MarkDuplicates的REMOVE_DUPLICATES=false,是行内默认配置。为什么标而不删?因为GATK HaplotypeCaller会在变异检测时基于概率模型处理duplicates信息,直接删除会在覆盖度低的区域引入额外的偏差。这些细节做管线设计时牵一发动全身,最好一开始就按社区最佳实践来。
3.4 第四步:BQSR碱基质量校正
BQSR(Base Quality Score Recalibration,碱基质量值校正)解决的是测序仪对碱基质量分数的系统性偏差问题。Illumina和DNBSEQ的base quality常常受测序循环数、dinucleotide上下文等因素影响,会导致质量分数与实际错误率不一致。
GATK BQSR分两步:
# 第一步:基于已知SNP位点计算校正模型 gatk BaseRecalibrator \ -R GRCh38.fa \ -I sample1.markdup.bam \ --known-sites dbsnp_146.hg38.vcf.gz \ --known-sites Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \ -O sample1.recal.table # 第二步:应用校正模型 gatk ApplyBQSR \ -R GRCh38.fa \ -I sample1.markdup.bam \ --bqsr-recal-file sample1.recal.table \ -O sample1.recal.bamBQSR对生殖系变异检测的最终结果提升不是"天翻地覆",但对低频突变、尤其等位基因频率在5%以下的位点,影响就体现出来了。肿瘤样本做低频变异检测时,BQSR这一步绝不能跳过。
4. 变异检测的核心逻辑:GATK HaplotypeCaller与多样本联合分析
处理完BAM之后就是变异检测。GATK HaplotypeCaller是当前全基因组SNV/InDel检测的社区标准,它比早期流程的统一基因分型(UnifiedGenotyper)在InDel附近的准确性有了质的提升。
4.1 HaplotypeCaller的算法逻辑
HaplotypeCaller不是简单地看单个碱基位置上有没有"不一致",而是采用**局部重新组装(local reassembly)**策略:在候选变异区域,把所有覆盖该区域的reads提取出来,做De Bruijn graph组装,构建出该区域可能的单倍型,再把每个read与这些单倍型进行比对,用概率方法判断哪个位点更可能是真实变异。
这个设计的好处是,在InDel附近,传统"逐位点比较"会被reads的比对歧义弄得很痛苦,而HaplotypeCaller直接绕过了这个问题——它能直接确定整个InDel的具体序列和边界。代价是计算量显著增加。
单样本变量检测命令:
gatk HaplotypeCaller \ -R GRCh38.fa \ -I sample1.recal.bam \ -O sample1.g.vcf.gz \ -ERC GVCF这里用-ERC GVCF模式,产出的是gVCF文件,记录每个位点的参考比对的基因型信息。若你只测了一个样本,直接用-O sample1.vcf.gz也行;但项目可能有后续新增样本,gVCF的产出可以方便做joint calling,强烈建议从一开始就按gVCF流程做。
4.2 Joint Calling与单样本Calling的差异
基因组测序项目中,如果同时有多个样本(比如一个家系、一个队列),最佳实践是分别对每个样本产出gVCF,然后统一用GenomicsDBImport合并,再做joint genotyping。这样做的好处:可以把所有样本一起对比参考序列,对低覆盖度位点的基因型判断更稳,也便于发现某些样本特异的低频信号。
# 建立genomicsdb gatk GenomicsDBImport \ -V sample1.g.vcf.gz -V sample2.g.vcf.gz -V sample3.g.vcf.gz \ --genomicsdb-workspace-path my_database \ -L intervals.list # 联合基因分型 gatk GenotypeGVCFs \ -R GRCh38.fa \ -V gendb://my_database \ -O cohort.vcf.gz需要注意-L intervals.list这个参数,GenomicsDBImport默认不接受没有interval的输入,你需要准备一个包含全部染色体的区间文件(比如从参考基因组的dict文件生成的chr1-22、chrX、chrY、chrM的列表),否则这步会直接报错。
单样本calling适合家系中先验样本数量少、不考虑横向比较的情况。但从管线可维护性的角度,我个人的选择是:只要数据量超过5个样本,就一律走gVCF+joint calling,后面无论追加样本还是下游新分析,都留有余地。
4.3 VQSR还是硬过滤
得到原始VCF之后,面对的问题是"怎么把假阳性过滤掉"。行内有两个主流方案:VQSR(Variant Quality Score Recalibration)和硬过滤(hard filter)。
VQSR的思路是利用高置信度的已知变异位点(如HapMap、1000 Genomes、dbSNP)作为训练集,用高斯混合模型对变异位点的多个特征(如QD、FS、MQ、MQRankSum等)打分,然后根据敏感度/特异性权衡取阈值。它需要至少30个样本来训练模型,否则容易出现过拟合。样本数不够时,走硬过滤更稳妥。
硬过滤常用参数(源于GATK best practice)大约是:
# SNV硬过滤 gatk VariantFiltration \ -V raw.vcf.gz \ -filter "QD < 2.0" --filter-name "QD_fail" \ -filter "FS > 60.0" --filter-name "FS_fail" \ -filter "MQ < 40.0" --filter-name "MQ_fail" \ -O filtered_snps.vcf.gz # INDEL硬过滤 gatk VariantFiltration \ -V raw_indel.vcf.gz \ -filter "QD < 2.0" --filter-name "QD_fail" \ -filter "FS > 200.0" --filter-name "FS_fail" \ -O filtered_indels.vcf.gz关于这个阈值有两点经验想说:QD<2.0比较保守,适合追求特异性的场景;如果你更看重敏感度(比如做罕见病研究,担心漏掉真变异),可以把QD阈值适当提高到2.5-3.0之间。FS(Fisher's strand bias)是衡量正负链支持是否均衡的指标,InDel的阈值比SNV高出一大截,就是因为InDel在链偏倚上天然比SNV更容易"看起来不平衡",直接照搬SNV阈值会把大量真实InDel过滤掉。
5. 变异注释与下游解读:从VCF到可用的生物学结论
到这里,变异已经拿到了,但一堆坐标和等位基因信息对生物学研究人员来说还不够——它们需要被翻译成人话和生物学意义。VCF中的位点需要注释到基因、转录本、氨基酸变化、人群频率数据库、致病性数据库等,这才是临床和科研真正依赖的信息。这一步也最容易看出一个生信人对"数据闭环"的把握能力。
5.1 常用注释工具:ANNOVAR与VEP
目前WGS变异注释的标配工具是ANNOVAR和Ensembl VEP。两者功能相似,但各有长短。
ANNOVAR的优势在于它的数据库体系非常成熟,注释结果简洁清晰,对大批量VCF的注释速度也快。用法大致是:
table_annovar.pl sample.vcf \ humandb/ \ -buildver hg38 \ -out sample.annovar \ -remove \ -protocol refGene,cytoBand,exac03,gnomad211_exome,gnomad211_genome,clinvar_20221231,avsnp150,dbnsfp41a \ -operation g,r,f,f,f,f,f,f \ -nastring . \ -vcfinput注意-protocol里包含了refGene(基因注释)、cytoBand(染色体区带)、gnomAD(人群频率)、ClinVar(临床意义)、dbSNP(常见位点)、dbNSFP(功能预测)等数据库。其中的-nastring .参数让注释不到的字段填充一个"."而不是留空,这个对下游筛选极其有用,因为留空字符和"."在Excel中的行为完全不同。
VEP的优势在于注释字段更标准,和Ensembl数据库更新同步,还支持自定义插件(如SpliceAI、CADD等)。如果你做的是面向临床解读的项目,建议VEP为主、ANNOVAR交叉验证。
5.2 注释字段里最容易被忽略的几个信息点
很多刚入门的人注释完就直接按"ExAC/gnomAD频率 < 0.01"筛变异,这个思路没问题,但会漏掉几个重要维度。
- CADD/PolyPhen/SIFT等功能预测:这些分数不是"真相",而是基于不同算法模型对变异的破坏性预测。它们容易误导人,因为它们对同一个位点的判断经常不一致,我的习惯是把这几个分数放一起看,而不是依赖单一分数。
- 剪接区域变异:WGS可以覆盖到内含子与剪接位点,这一部分对临床意义极大。约15%的致病性点突变发生在剪接供体/受体位点——这些位点在WES上往往覆盖很差,而在WGS中可以稳定检测到。
- 非编码区的调控注释:如ENCODE的候选顺式调控元件(cCREs)、Roadmap的表观遗传标记,这些对解释非编码区变异非常关键。目前对非编码区变异的解读仍处于早期,但使用RegulomeDB或者OpenTarget Genetics等工具查找还是很有帮助。
5.3 结构变异与拷贝数变异的检测思路
WGS的结构变异检测是它的核心优势。工具选择上,常用的有Manta(侧重InDel和SV的breakpoint精确定位)、DELLY(支持deletion、duplication、inversion、translocation四类经典SV)、Lumpy(利用split-read和read-depth进行综合检测)。实际操作中,大家通常会合并多个工具的预测结果,而Manta是这里面灵敏度与精度平衡较好的标准选择。
CNV分析通常是利用覆盖度(read depth)推断拷贝数,常用工具是CNVkit、Sequenza或PureCN。对于肿瘤WGS样本,结合BAM的等位基因频率和杂合性丢失(LOH)信息,可以推导更复杂的allele-specific拷贝数,PureCN和ASCAT是其中的代表。对于生殖系样本,如果有正常对照配对,则CNVkit更适合研究。
6. 肿瘤WGS的特殊处理:纯度、倍性与克隆结构
肿瘤样本的WGS分析,分析难度完全不是一个量级。肿瘤组织里混着正常细胞,肿瘤细胞内部还有亚克隆(subclones),每个亚克隆的拷贝数和变异谱又不一样。如果不做纯度/倍性估计,后续变异检测、CNV推断、克隆结构分析都会失真。
6.1 纯度与倍性估计对检测的影响
肿瘤样本中,如果肿瘤细胞占比只有40%——也就是说变异等位基因频率(VAF)理论上限是0.2(纯合变异)或0.1(杂合变异)——此时直接套用生殖系的变异检测阈值就会漏掉大量真正的体细胞突变。所以第一步是估计肿瘤纯度(purity)和倍性(ploidy)。
工具层面,Sequenza通过分析BAM的B-allele frequency(BAF)和read depth的比值来估计纯度和倍性,输出肿瘤基因组每个染色体臂的拷贝数状态;ASCAT同样是基于BAF和LRR(log R ratio)做估计。这类工具输出的绝对拷贝数(total copy number与major/minor allele count)比单纯的相对拷贝数更有价值,你可以直接从结果中判断某个片段是杂合性缺失(LOH)还是等位基因不平衡。
6.2 体细胞变异检测的配对分析
肿瘤WGS的体细胞变异检测必须要有配对正常样本(通常是外周血或癌旁组织),否则无法区分"种系固有变异"和"肿瘤特有变异"。Musical架构下,Mutect2是最常用的体细胞SNV/InDel检测工具,它通过tumor-normal配对模型,将正常样本中存在的变异从肿瘤样本中扣除,同时利用panel-of-normals(PoN)校正批次性测序错误。
gatk Mutect2 \ -R GRCh38.fa \ -I tumor.bam \ -I normal.bam \ -tumor TUMOR_SAMPLE_NAME \ -normal NORMAL_SAMPLE_NAME \ --germline-resource af-only-gnomad.hg38.vcf.gz \ --panel-of-normals pon.vcf.gz \ -O tumor_normal_m2.vcf.gz--panel-of-normals和--germline-resource这两个输入文件对过滤测序平台系统错误意义重大。PoN能显著压低"所有样本都会出现的假阳性位点"——这类假阳性往往与探针/平台偏好性相关,纯靠tumor-normal配对根本去不掉。没有PoN时,Mutect2的假阳性率会明显上升。
6.3 肿瘤克隆结构与进化推断
体细胞变异做出来之后,很多项目要回答的问题是:这些突变是同一个克隆携带的,还是肿瘤内部存在不同亚克隆?哪个克隆是祖先?这涉及肿瘤克隆进化推演。常用工具是PyClone或SciClone,它们的核心逻辑是基于每个位点的VAF以及该位点对应的拷贝数状态,把共变的位点聚类成不同的克隆群,再进一步推断克隆间的出现顺序。
这一块的分析内容比较深,这里只补充一个关键点:克隆推断依赖的VAF值,必须经过纯度与拷贝数校正(即从观察VAF换算为"肿瘤细胞中携带该突变的细胞占比"),否则克隆聚类会被系统性高估或低估。
7. WGS数据量、存储与算力规划:做生信之前先做"基建"
坦白讲,很多人流程都跑通了,最后被"数据存不下"或"服务器算不动"卡住。WGS的数据体量非常实在:一个样本的fastq约60-80GB,BAM约80-120GB,gVCF约1-2GB,VCF约100MB左右。一个20个样本的WGS队列,从原始数据到最终分析结果,中间产物加总通常在2TB上下,如果没有合理的存储策略,很容易出问题。
7.1 存储规划建议
- fastq是原始数据,建议双副本存档(本地+冷备/云对象存储),因为重新测序的成本远高于存储成本。
- BAM是分析中间产物,建议保留至少一份。如果你的分析流程需要反复调整,保留BAM能避免重新比对(那是最大计算量环节)。
- gVCF/VCF是分析产出,体积小又是下游所有分析的基础,必须长期保存。
- 中间临时文件(如未排序的BAM、中间分割文件)可以及时删除,不用心疼。
7.2 算力规划与调优
全基因组比对、排序和去重这三个步骤是算力消耗的大头。一个30x WGS样本,如果用32核CPU的机器跑典型GATK流程,从fastq到变异检测完成,大约需要6-10小时。批次处理多份样本时,建议优先考虑云上批量计算服务,而不是本地单机硬扛,因为WGS流程天然可并行化——不同样本可以同时跑比对,样本内部不同染色体可以并行处理。
线程参数也有讲究。BWA-MEM的-t不是越大越好,它受限于磁盘I/O和内存带宽;samtools sort的-m参数(每个线程的最大内存)要结合机器总内存合理设置,否则容易OOM。对于标配64核、256GB内存的机器,我常用的配置是bwa mem -t 32,samtools sort -@ 16 -m 2G。
7.3 容器化与流程管理
当前流程管理的行业事实标准是Snakemake与Nextflow。两者都支持容器化执行(Docker/Singularity),能保证流程在不同机器上的可复现性。对于WGS,我建议使用nf-core/sarek(Nextflow生态),它整合了从BAM预处理到Variant calling再到注释的完整WGS/WES流程,自带多个标准工具与最佳实践参数配置。比从零搭管线省力得多。
示例:用Sarek跑WGS germline流程
nextflow run nf-core/sarek \ -profile docker \ --input samplesheet.csv \ --genome GATK.GRCh38 \ --tools haplotypecaller,strelka,manta \ --outdir results \ -resume它在底层已经把BQSR、gVCF、joint calling等步骤串起来了,速度、资源管理、错误重试都做了优化。直接用这套流程比自己手搓更稳。
8. 实战中的几个"坑位":覆盖度不均、性别染色体与重复区域
最后这部分没有系统的逻辑链条,纯粹是这些年踩过的坑,每条都对应过真实的返工和加班。
8.1 覆盖度不均:做得深不如覆盖得均匀
很多人在设计WGS项目时只关心"平均测序深度达到30x",但平均深度是一个被严重高估的指标。基因组不同区域的比对效率差异巨大——GC含量极端区域(尤其GC>70%或GC<30%)、着丝粒区域、端粒区域、rDNA重复簇等,测序仪和比对工具在这里都容易翻车。真正的判断指标是"基因组中达到20x覆盖度的区域占比",行业经验值是超过90%(也就是常说的'>90% of bases at 20x')。使用bedtools或mosdepth可以方便地统计这个值:
mosdepth --by 500 sample sample.recal.bam # 输出sample.per-base.bed.gz,用awk/python统计>=20x的Base比例以我的经验,一份看起来"平均30x"的数据,如果超过8%-10%的区域连15x都不到,那这份数据做低频变异检测会心虚——尤其当你关注的位点恰好落在低覆盖区域时。
8.2 性染色体变异的处理不能想当然
这是一个特别容易出问题又特别容易被忽视的点。女性和男性样本的chrX、chrY覆盖度天然不同,直接套用常染色体的过滤阈值会导致一堆假阳性/假阴性。比如chrY在女性样本中原则上比对不到(除非参考基因组有同源区域),但实际比对会有零星的Mappable reads mapped到chrY的假阳性;男性样本的chrX是半合子,杂合性判断与常染色体完全不同。建议在变异检测与过滤阶段就按染色体类型(autosome/chrX/chrY/mtDNA)分组考虑阈值,不要一刀切。
8.3 高重复区域的假阳性:SETD2、C4等例子
人类基因组中富含节段重复(segmental duplication)区域,如HLA区域、C4补体基因簇、PMS2基因的假基因区域等。这类区域在比对时会出现"多个位置都能match"的情况,BWA-MEM会给出MAPQ较低的比对结果。如果硬是把这些区域纳入分析,并且不加过滤,产出的"变异"有很大概率是比对噪声。
处理策略包括但不限于:分析前对低MAPQ(比如MAPQ<20或<30)的reads进行过滤;对高重复基因(如PMS2、CHEK2的部分外显子)使用专门的repeat-aware分析工具或已知的难处理区域列表(如UCSC可mappability track)来标注结果。很多人第一次拿到PMS2基因的"致病突变"时兴奋了一下,然后发现它在假基因上的同源性比真实区域还高。
8.4 参考基因组版本:统一才是硬道理
hg19还是hg38?这个问题如果项目启动时没有约定好,后面会引发连锁灾难。参考基因组的坐标体系不一样,变异注释和临床解读时与公共数据库的比对就会错位。早期很多数据库以hg19为主,但现在gnomAD、ClinVar等主流数据库都已全面支持hg38,且hg38对许多区域的填补和修正做得更好。我的建议是:新项目一律用hg38(GRCh38),如无特殊理由不用hg19。历史遗留的hg19数据,建议用CrossMap或liftOver做坐标转换,但转换后的数据一定要做抽样检查,因为liftOver在某些复杂区域(尤其InDel附近)会丢位点。
9. 写在一份WGS报告末尾的话
之前带过的一个学生问我:"WGS数据分析流程学完,是不是就能做任何组学分析了?"我当时的回答是:工具会用了,但离'会做'还差一个维度——你得知道你手上的数据是怎么来的、哪些环节可能在撒谎、哪些区域要打折扣看待。测序数据本质上是一个概率推断问题,每一个结论背后都建立在若干假设之上;真正从事这个行业越久,越觉得生信的核心不是命令敲得有多快,而是能不能在数据与生物学之间找到一个可靠的桥梁。
再分享一个习惯:无论跑什么样本,处理完之后我都会挑几个已知阳性的位点(比如你们项目内部验证过的突变)作为质控基准,确认流程能稳定检测到它们,再批量推进。这个习惯帮我拦下过至少三次由上游样本标签错乱引发的灾难级事故。
WGS这条路不短,从实验台到服务器,从FASTQ到解读报告,每个环节都有坑。这篇内容把我能想到的、实操中最关键的信息都写进来了,希望对正在做或者准备做WGS的你有些帮助。你也可以先把这篇收藏起来,等实际跑到某一步卡住时再回来翻对应的小节——很多细节不是在读的时候记住的,是在踩坑的时候才真正理解的。