☰
全基因组测序(WGS)分析流程详解:从FASTQ到VCF的实战指南
2026/10/5 13:32:42 网站建设 项目流程

WGS这个名头,做生信的朋友应该都不陌生,但真正要把全基因组测序的数据从头到尾跑一遍,从FASTQ一路处理到VCF、再到注释结果,中间的门道远不止“软件挨个跑一遍”那么简单。我自己从最开始拿公共数据练手,到后面接手真实样本的WGS分析任务,踩过的坑、推翻重来的流程都不少。这篇文章就把我实际跑WGS项目时的完整思路、工具选型逻辑、每一步的具体操作和参数调整经验整理出来,希望对刚入坑生信、或者准备从WES转WGS的朋友有点帮助。

1. WGS到底是什么,为什么值得做

1.1 WGS与WES/RNA-seq的定位差异

全基因组测序(Whole Genome Sequencing,WGS)顾名思义,是对一个物种的整个基因组进行测序,覆盖范围包括所有编码区和非编码区。相比之下,全外显子测序(WES)只捕获约1%-2%的蛋白编码区域,RNA-seq则主要关注转录组层面的表达情况。

这三者最容易让人混淆的就是WGS和WES。一句话总结:如果想知道“基因序列里到底哪里有变异”,WGS是最彻底的方式,因为它不分区域、不做捕获富集,把基因组里能测的都测了。外显子组测序虽然便宜,但测序原理决定了它必须先经过探针杂交捕获,这个过程会引入捕获偏好性,一些GC含量极端区域的覆盖度往往不理想。

举个实际例子:一个肿瘤样本WES只覆盖了约30Mb的区域,而WGS覆盖约3Gb,两者数据量相差近百倍。分析WGS时,BAM文件大小、VCF文件行数、计算资源占用都要按数量级去重新估算。

1.2 选择WGS时的成本和数据量考量

成本是绕不开的话题。WGS的测序单价这些年降得很快,但即便按个人全基因组测序市场价3000-5000元来计算,一个样本30X深度的原始数据就有大约90-100Gb的FASTQ文件。如果做的是肿瘤配对样本(肿瘤+对照),轻松就是200Gb的原始数据起步。

这里给一个数据量的估算公式,方便大家做项目规划:

  • 人类基因组大小约为3.2Gb(含参考基因组的所有染色体)
  • 30X深度的总测序量约为 3.2Gb × 30 ≈ 96Gb
  • 测序仪下机数据转换成FASTQ,去除低质量reads和接头后,实际有效数据约占总测序量的85%-95%
  • BAM文件大小约为FASTQ的1.5倍左右,排序去重后再压缩,通常和FASTQ大小相当

熟悉这个量级后,就知道为什么WGS项目必须用高性能计算集群或者高配置工作站来跑,而不是指望普通笔记本能撑住。存储空间也是大头,中间文件(BAM、GVCF)千万别随便删,后面统计分析阶段很可能还要回溯。

1.3 WGS的核心应用场景

从实际应用来看,WGS在以下几个场景里是刚需:

  • 罕见病诊断:不仅能看到SNP和InDel,还能检测拷贝数变异(CNV)、结构变异(SV),对于临床上找不到原因的疑似遗传病,WGS的检出率明显高于WES
  • 肿瘤体细胞变异检测:通过肿瘤样本与正常对照配对测序,可以全面刻画点突变、插入缺失、拷贝数改变以及染色体结构重排,为靶向治疗和免疫治疗提供依据
  • 微生物基因组研究:细菌、病毒全基因组测序用于分型、耐药基因鉴定、暴发溯源
  • 群体遗传与进化分析:大批量样本的低深度WGS可以在控制成本的同时获得全基因组范围的变异信息

我自己的项目背景是肿瘤WGS,所以后面讲流程时会更侧重这个方向的细节。不过整体分析框架在人类WGS里是通用的。

2. WGS分析流程设计与工具选型拆解

2.1 标准分析pipeline总览

一个典型的WGS生殖系变异分析流程可以用下面这张表来概括:

步骤主要工具输入输出核心目的
数据质控FastQC, MultiQC, fastpFASTQ干净FASTQ、质控报告去除低质量碱基和接头
序列比对BWA-MEM2 / BWA-MEM干净FASTQSAM/BAM将reads定位到参考基因组
排序与去重Samtools, Picard MarkDuplicatesBAM排序去重后的BAM消除PCR重复对变异频率的影响
碱基质量校正GATK BaseRecalibrator / ApplyBQSRBAM校正后的BAM校正测序仪系统误差
变异检测GATK HaplotypeCaller / DeepVariantBAMVCF / GVCF检出SNP和InDel
变异过滤GATK VariantRecalibrator / 硬过滤VCF过滤后的VCF去除假阳性变异
注释与解读ANNOVAR / SnpEff / VEPVCF注释后的表格给变异加上基因、功能影响等注释

这里要注意:这个流程是“生殖系变异”的标准分析路径。如果是肿瘤体细胞变异,第二步到第五步仍然类似,但变异检测环节会更推荐用Mutect2这类工具做配对分析,过滤逻辑也完全不同。

2.2 比对工具为什么选BWA

BWA-MEM几乎成了WGS比对的“默认选项”,原因是它兼顾了速度和准确性,在人类基因组这种复杂重复结构背景下依然稳定。BWA-MEM2进一步对算法做了优化,利用SIMD指令集加速,实测下来比原版快大约2倍,内存占用也差不多。

可能有人会问:为什么不用Bowtie2或者STAR?Bowtie2在转录组比对里用得比较多,但它在处理长读长(150bp PE)和indel时不如BWA-MEM;STAR虽然非常快,但它本质上是为RNA-seq设计的剪接比对工具,在WGS这种DNA测序场景里不是最优选择。工具选型永远是:先明确场景,再对比工具。

比对时还值得注意参考基因组的选择。人类WGS分析现在主推GRCh38,老旧的GRCh37/hg19基因注释和坐标系统差异很大,后续如果要做临床注释,尽量从源头就用GRCh38,避免转换时引入坐标错误。实际工作中我见过不少项目还在用hg19,倒不是说完全不行,但最后在变异注释和跨项目比较时会很痛苦。

2.3 GATK HaplotypeCaller与DeepVariant的对比

变异检测是这个pipeline里最核心、也最能体现分析者水平的一步。

GATK HaplotypeCaller是一个经典的“组装型”变异检测器,它会对每个区间从头组装局部单倍体,再用PairHMM算法计算每个单倍体的似然值,最后输出最可能的基因型。优点是成熟、稳定、社区资料多,在Illumina平台、WGS场景下有大量验证数据支撑。

DeepVariant则是Google基于深度学习的方法,把“变异检测”变成“图像识别”问题——将比对结果编码成带有颜色的图像,再用神经元网络做分类。实测下来,DeepVariant在WGS数据上的精确率和召回率通常优于GATK流程,尤其在InDel检测上提升明显。

但我个人的建议是:如果想跑通一套流程、方便在社区里寻找帮助,优先用GATK流程。想追求极限准确性、算力充裕(尤其是GPU资源),可以拿DeepVariant做交叉验证。两个工具结果取交集和并集,有时能发现一些单工具检不出的边缘变异。

2.4 注释工具怎么选

变异注释是连接生信分析和临床解读的桥梁。

  • ANNOVAR:中文作者开发,学术免费,功能全面,输出形式灵活,可以做基因注释、区域注释、数据库注释(dbSNP、1000G、ExAC等)。缺点是需要自己下载和维护数据库文件
  • SnpEff:Java写的,速度快,内置数据库丰富,注释结果自带变异影响等级,适合通量较大的项目
  • VEP(Variant Effect Predictor):Ensembl官方出品,注释信息详尽,支持插件扩展,但运行速度相对慢一些

我的倾向是:如果只是快速批量注释,SnpEff就够了;如果想给临床报告或者科研论文提供更丰富的注释信息,推荐用ANNOVAR或者VEP。很多人习惯把SnpEff和ANNOVAR都跑一遍做交叉验证,注释结果一致性越高,后续分析越放心。

3. 实操记录:从一个30X人类WGS样本开始

3.1 环境准备与参考基因组配置

我自己用的跑WGS的环境是Linux服务器,内存64GB以上,核心数16个以上,这样单样本跑起来不会太吃力。软件方面优先使用conda环境,方便管理依赖版本:

# 创建conda环境并安装核心工具 conda create -n wgs python=3.9 conda activate wgs conda install -c bioconda fastqc fastp multiqc bwa samtools picard gatk4

参考基因组这里强调一下:千万不要直接在网上下载一个hg38.fa就开始比对,标准的做法是在GATK官方提供的资源包基础上准备基因组字典文件、索引文件。

# 下载GATK bundle资源中的GRCh38(通常包含primary assembly和decoy序列) # 例如从官方FTP或GCP bucket获取 Homo_sapiens_assembly38.fasta # 生成BWA索引 bwa index Homo_sapiens_assembly38.fasta # 生成.fai索引 samtools faidx Homo_sapiens_assembly38.fasta # 生成dict字典文件(GATK需要) gatk CreateSequenceDictionary -R Homo_sapiens_assembly38.fasta

准备工作大概需要一两个小时,取决于磁盘IO和CPU性能。索引文件会占据不少磁盘空间,BWA的索引加.fai、dict,整体可能多出几个Gb的占用,做好心理准备。

3.2 原始数据质控:别上来就比对

拿到下机数据之后,我习惯先跑一个FastQC和MultiQC的质性检查。多线程条件下,单样本大概10分钟左右就能看到结果。

# 创建原始FASTQ和输出目录 mkdir -p raw_data fastqc_result clean_data # 查看原始数据量 ls -lh raw_data/ # 通常会看到两个文件:sample_R1.fastq.gz 和 sample_R2.fastq.gz # FastQC检查原始数据 fastqc -t 8 -o fastqc_result raw_data/sample_R1.fastq.gz raw_data/sample_R2.fastq.gz # 聚合QC报告 multiqc -o fastqc_result fastqc_result

FastQC报告里我主要关注这几项:Per base sequence quality(Q30是否达到80%以上)、Per sequence GC content(有没有异常双峰)、Overrepresented sequences(是否有接头残留或纯净重复)、Adapter Content(接头率是否高)。

如果接头率偏高或者质量分数低,下一步就要做修剪。fastp是我目前用下来最快最顺手的工具,一步到位做质量过滤、接头切除、低复杂度序列过滤,而且还会生成一个很直观的HTML报告。

fastp \ -i raw_data/sample_R1.fastq.gz \ -I raw_data/sample_R2.fastq.gz \ -o clean_data/sample_R1.clean.fastq.gz \ -O clean_data/sample_R2.clean.fastq.gz \ --detect_adapter_for_pe \ --thread 16 \ --html clean_data/sample_fastp.html

不少教程会同时保留Trimmomatic做对比,实际上在新手项目里我建议就选fastp,配置文件简单,参数透明,跑完能直接看到详细统计。如果测序质量本身已经很理想,也可以跳过修剪直接比对,但为了后续结果纯净,我一般还是会过一遍fastp。

3.3 比对、排序、去重的完整过程

这一步是整个流程中最吃CPU和内存的环节,也是我一开始最容易“卡壳”的地方。为什么这么说?因为当你把BWA-MEM的输出直接重定向成BAM文件时,会发现最终得到的文件是未排序的,后续步骤全都没法直接用。我最初跑的时候不懂这个,直接在比对命令后面接了一个标记重复的步骤,结果报错报了半小时。

正确做法是:比对、按坐标排序、标记重复,分步进行。

# 设置参考基因组路径 REF=Homo_sapiens_assembly38.fasta SAMPLE=sample # 步骤1: BWA-MEM2比对(推荐MEM2,速度更快) bwa-mem2 mem \ -t 16 \ -R "@RG\tID:${SAMPLE}\tSM:${SAMPLE}\tPL:ILLUMINA\tLB:${SAMPLE}" \ ${REF} \ clean_data/${SAMPLE}_R1.clean.fastq.gz \ clean_data/${SAMPLE}_R2.clean.fastq.gz \ > ${SAMPLE}.sam # 步骤2: SAM转BAM并排序 samtools view -bS -@ 16 ${SAMPLE}.sam | \ samtools sort -@ 16 -m 4G -o ${SAMPLE}.sorted.bam # 步骤3: 标记PCR重复 gatk MarkDuplicates \ -I ${SAMPLE}.sorted.bam \ -O ${SAMPLE}.dedup.bam \ -M ${SAMPLE}.dedup.metrics.txt \ --CREATE_INDEX true # 步骤4: 生成比对统计 samtools flagstat ${SAMPLE}.dedup.bam > ${SAMPLE}.flagstat.txt

为什么必须加上-R参数?因为它定义了Read Group信息,这是GATK后续做Base Quality Score Recalibration(BQSR)和变异检测过程中区分样本来源的关键。如果多个样本合并分析时没有Read Group信息,后面会乱成一锅粥。

MarkDuplicates这个步骤初看起来“只”是为了去掉PCR重复,但在WGS数据里它直接决定了覆盖深度的准确性。一个30X的WGS样本,PCR重复率一般在5%-15%之间,重复率过高说明文库复杂度不够,后续的变异频率计算就会失真。

另外,关于-m 4G这个参数的设置也要提一句:它控制的是samtools sort时每个线程的最大内存,不是总内存。如果机器内存大、线程多,可以适当提高这个值来加速排序;反之,如果总内存只有32GB、线程数16,那么-m 2G更稳妥,否则容易直接内存溢出。

3.4 BQSR碱基质量校正:不可跳过的严谨步骤

比对去重完成后,建议执行BQSR(Base Quality Score Recalibration)。这一步的目的是修正测序仪在不同基序、不同位置上的系统误差。在WGS这种追求极致灵敏度的数据里,它能让变异检测的假阳性明显下降。

# 已知位点数据库(用GATK资源包自带的dbsnp和1000g gold standard) DBSNP=resources_broad_hg38_v0_Homo_sapiens_assembly38.dbsnp138.vcf GOLD_1000G=resources_broad_hg38_v0_1000G_phase1.snps.high_confidence.hg38.vcf.gz # 步骤1: 生成线性模型 gatk BaseRecalibrator \ -R ${REF} \ -I ${SAMPLE}.dedup.bam \ --known-sites ${DBSNP} \ --known-sites ${GOLD_1000G} \ -O ${SAMPLE}.recal.table # 步骤2: 应用校正到BAM文件 gatk ApplyBQSR \ -R ${REF} \ -I ${SAMPLE}.dedup.bam \ --bqsr-recal-file ${SAMPLE}.recal.table \ -O ${SAMPLE}.recal.bam

很多人觉得BQSR跑起来慢,耗时和变异检测差不多,干脆就跳过了。但在低深度WGS或者肿瘤样本中,跳过这一步可能让一些位点因碱基质量的系统性偏差被错误检出。我个人建议:标准分析不要省这一步骤。如果要优化速度,至少保证比对和去重是正确的。

3.5 变异检测:HaplotypeCaller实战

变异检测这一步,GATK HaplotypeCaller可以按两种模式运行:

  • 传统模式:直接输出VCF,适合单样本分析
  • GVCF模式:输出每个位点的基因型信息(包含位点调用置信度),适合多样本联合分析时先各出各的GVCF,再用GenomicsDBImport合并

实际项目中,即使只有一个样本,我也更推荐输出GVCF,因为后续可以随时和其他样本合并,不需要重新跑比对和变异检测。

# GVCF模式 gatk HaplotypeCaller \ -R ${REF} \ -I ${SAMPLE}.recal.bam \ -O ${SAMPLE}.g.vcf.gz \ -ERC GVCF \ --native-pair-hmm-threads 16

注意:HaplotypeCaller非常吃CPU和内存,单样本30X WGS在这个步骤大概需要6-10个小时,线程给足、内存足够的情况下能明显提速。在此期间尽量避免跑其他重负载任务,否则容易因为资源竞争导致运行极慢。

GVCF文件生成后,如果是单个样本,直接用GenotypeGVCFs获得最终VCF:

gatk GenotypeGVCFs \ -R ${REF} \ -V ${SAMPLE}.g.vcf.gz \ -O ${SAMPLE}.vcf.gz

如果是多个样本,需要先用GenomicsDBImport批量导入GVCF再做联合基因分型。这一步数据库的磁盘占用也很大,建议预留200Gb以上的空间。

3.6 变异过滤:VQSR还是硬过滤

变异过滤是决定最终结果可信度的关键环节。GATK里主要有两种过滤方式:

VQSR(Variant Quality Score Recalibration,变异质量分数校正):利用已知的高置信度数据集训练模型,用机器学习的思想评估变异质量。它的前提是需要足够多的变异位点(通常要有8000个以上SNP位点),因此最适合WGS,WES因为位点数不够、模型容易出错。

硬过滤(Hard filter):按QD、FS、MQ、MQRankSum等指标直接设一个阈值过滤。在样本少、位点数不足的情况下更稳健。

标准WGS分析里,如果样本量不够做VQSR,我建议先做硬过滤:

# 分离SNP和InDel,分别做硬过滤 gatk SelectVariants -R ${REF} -V ${SAMPLE}.vcf.gz --select-type-to-include SNP -O ${SAMPLE}.snp.vcf.gz gatk SelectVariants -R ${REF} -V ${SAMPLE}.vcf.gz --select-type-to-include INDEL -O ${SAMPLE}.indel.vcf.gz # SNP硬过滤(GATK推荐阈值) gatk VariantFiltration \ -R ${REF} \ -V ${SAMPLE}.snp.vcf.gz \ --filter-expression "QD < 2.0 || FS > 60.0 || MQ < 40.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" \ --filter-name "SNP_HARD_FILTER" \ -O ${SAMPLE}.snp.hardfilter.vcf.gz # InDel硬过滤 gatk VariantFiltration \ -R ${REF} \ -V ${SAMPLE}.indel.vcf.gz \ --filter-expression "QD < 2.0 || FS > 200.0 || ReadPosRankSum < -20.0" \ --filter-name "INDEL_HARD_FILTER" \ -O ${SAMPLE}.indel.hardfilter.vcf.gz # 合并回到一个VCF gatk MergeVcfs -I ${SAMPLE}.snp.hardfilter.vcf.gz -I ${SAMPLE}.indel.hardfilter.vcf.gz -O ${SAMPLE}.filtered.vcf.gz

硬过滤之后的VCF里,每个位点的FILTER列会标记是PASS还是被过滤掉了,后续分析通常只保留PASS的位点。

如果样本量足够(比如几百例WGS),VQSR的过滤效果会更平滑,能保留更多低但可信的边缘位点,不过计算量也不小。这个可以根据项目和计算资源自行判断。

3.7 变异注释:从位点到生物学意义

前面得到的是一个包含数百万行变异信息的VCF,但纯粹看坐标、等位基因,仍然不知道它有什么生物学意义。注释这一步就是给变异“贴标签”。

以ANNOVAR为例,一个最常见的注释流程如下:

# 下载或准备数据库(选核心的几个:refGene, dbnsfp, gnomad, clinvar) # 假设数据库目录为 annovar_db/ # 将VCF转换为ANNOVAR格式 annovar/convert2annovar.pl \ -format vcf4 \ ${SAMPLE}.filtered.vcf.gz \ > ${SAMPLE}.avinput # 做基因注释和数据库注释 annovar/annotate_variation.pl \ -geneanno \ -dbtype refGene \ -buildver hg38 \ ${SAMPLE}.avinput \ annovar_db/ annovar/annotate_variation.pl \ -filter \ -dbtype clinvar_20221231 \ -buildver hg38 \ ${SAMPLE}.avinput \ annovar_db/

注释结果会生成一个.exonic_variant_function文件,记录了每个变异落在哪个基因的外显子区、是错义突变还是同义突变,还会给出ClinVar里是否有已知致病性记录。这些信息对后续筛选候选致病位点非常有帮助。

如果是肿瘤样本,还建议加上COSMIC数据库以及肿瘤相关的热点变异注释(OncoKB等),这样可以直接筛选驱动突变。

4. 常见问题与排查技巧实录

4.1 比对率异常低?先查这四件事

WGS分析最让人“心态爆炸”的,就是跑到比对这一步,发现比对率不到80%,甚至只有60%。这种情况首先不要急着怀疑测序样本,按下面的顺序排查:

第一,参考基因组版本与reads是否匹配。如果生物样本是来自人类,但参考基因组错下了细菌的序列,比对率必然惨不忍睹。确认一下参考基因组源头和物种信息。

第二,FASTQ是否被污染或链路错误。比对之前,用FastQC看一下是否有明显的接头残留、未知序列比例过高。有时候样本标签写反了,R1和R2搞混,也会严重影响比对率。

第三,是否有过多无参考序列的reads。一些公共数据里的reads可能来自人工构建序列、载体或线粒体,这些序列如果在参考基因组里没被收录,会被判为未比对。WGS正常比对率在95%以上,98%左右是健康状态,如果低于90%,问题就比较大。

第四,覆盖度和深度是否低估。有时候比对率看似正常,但某个染色体的覆盖深度特别低,这时候要考虑是不是样本存在染色体异倍体或者参考基因组构建质量问题。

我记得有一个项目,样本来自有一定血缘关系近亲的个体,比对率本身很正常,但InDel的假阳性率特别高,后来排查发现是参考基因组的某些区域存在较多的同源序列,导致BWA-MEM2的map质量偏低。这时候可以适当调高比对质量阈值,比如过滤MAPQ < 20的reads,能降低一部分假阳性。

4.2 重复率高得离谱,问题可能在文库片段

PCR重复率高,有时候不见得是分析流程错了,而是文库构建出了状况。WGS文库构建时,如果起始DNA量过低,PCR循环数太多,或者片段化长度过短,都容易造成重复率飙升。正常WGS文库的重复率在5%-15%之间,如果超过20%,就有必要重新评估文库质量了。

在分析侧,如果你的重复标记阶段用samtools markdup,功能和Picard MarkDuplicates很接近,但某些情况下对双端reads的确认逻辑还是略有差异。为了和上游数据分析保持一致性,通常推荐直接使用GATK中的MarkDuplicates,毕竟它和HaplotypeCaller在同一套官方流程里。

这里也给出一个排查重复率的快捷命令:

# 查看重复率指标 grep -E "Unknown Instrument|PCT_DUPLICATION" ${SAMPLE}.dedup.metrics.txt

PCT_DUPLICATION一般以小数表示,比如0.10就是10%的重复率,超过0.25就要高度警惕。

4.3 变异检测结果“忽多忽少”的常见原因

同样的样本、同样的数据,换了一个流程版本,变异数量可能有几百上千的差异。这里最容易出问题的点:

第一,GATK版本差异。GATK3和GATK4的算法虽然同源,但参数默认值、模型版本都有差别。GATK4里HaplotypeCaller默认不再做物理相位和并行变异检测的同时处理,OpenMP线程参数也不一样,所以结果会有细微差异。跨项目比较时尽量锁定GATK4.x的某个具体版本。

第二,GVCF合并方式的差异。多样本分析时,如果用传统CombineGVCFs和GenomicsDBImport分别做合并,结果也可能略有出入。官方推荐GenomicsDBImport,不仅速度快,而且更适合大数据集。

第三,参考基因组的“decoy”序列。GRCh38的primary assembly和加入decoy(病毒、病毒载体等)序列的版本,比对时多数reads会落在真正的染色体上,但decoy序列能吸附一些无法比对的reads,从而避免它们被错误分配到其他地方。如果和之前的项目参考不一致,结果对比时会有大量假阳性差异。

我自己在实操中的经验是:做WGS分析之前,把参考基因组、软件版本、GATK bundle数据库版本全部固定下来并写进README,这是标准化的第一步。

4.4 计算资源不够时的替代方案

不是所有人都有高性能计算集群。如果只有一台普通工作站(比如16核、64GB内存),跑30X人类WGS依然可行,但要对流程做几个策略性调整:

  • 比对和排序之间用管道连接,减少中间SAM文件占用的磁盘空间
  • MarkDuplicates阶段使用--CREATE_INDEX true减少一步额外建索引的时间
  • HaplotypeCaller指定--native-pair-hmm-threads为4-8,避免线程数过大导致内存爆炸
  • 在磁盘空间紧张的情况下,可以考虑在比对完成后删除原始FASTQ(如果存储了下机数据的备份),保留BAM和VCF就足够大多数下游分析使用了

如果连这一步都吃力,还有一个思路是做低深度WGS(5-10X),虽然每个位点的检测敏感性有所下降,但在群体层面做频率分析和CNV,要比跑到一个高深度样本却卡死在计算资源上靠谱得多。这个方案在遗传学入门项目里特别常见。

4.5 一些容易被忽略的细节问题

多平台数据合并:有的样本在Illumina不同型号的测序仪上跑了多个lane,合并FASTQ时要确保两个lane的Read Group信息保持区分,不能混为一谈。合并之后最好检查一下每个RG的reads数量。

chr前缀不一致:VCF注释过程中,有时会遇到参考基因组染色体命名不一致的情况——比如有些版本用1,有些用chr1。如果不定好统一标准,注释结果里可能有一半位点匹配不上。我的习惯是全程使用Ensembl/GATK资源包的命名(如1、2、3、X、Y),一旦用chr1,后面的GATK数据库也要配套。

重复样本的批次效应:如果项目里不同批次的结果差异明显,建议先做一个PCA分析(基于VCF的常见变异),看看样本是否按照批次聚类。如果是,需要在后续关联分析中把批次作为协变量,或者重新检查文库构建和测序条件。

写在最后的一点经验

WGS分析跑通了之后,我最大的体会是:流程本身其实并不复杂,难的是在每个环节都保持对数据质量的敏感。比对率、重复率、插入片段大小分布、VCF的Ti/Tv比值(转换/颠换比,人类WGS通常在2.0以上)……这些指标每一处都透露着数据质量的信号,只看你愿不愿意多花几分钟去检查。

最开始跑WGS时我也犯过“一键跑完脚本,然后拿着VCF直接做下游”的错误,结果到注释阶段才发现比对环节里混入了大量适配片段,只能从头再来。从那以后我养成了一个习惯:每一大步完成之后,先看一眼统计指标再往下走。虽然看起来效率低了一点,但回头排查问题节省的时间,远比这些零碎检查多得多。

如果你正准备开始自己的第一个WGS项目,先把参考基因组和软件版本固定好,再照着完整流程跑一遍,仔细看看每一步的产物和统计指标。别急着追求速度,把每个环节的原理搞清楚,后续做任何变异筛选、比较分析都会顺手很多。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询