生信分析初测实战：从FASTQ质控到变异检测的完整流程

简介：近年来生物信息岗位面试愈发注重Linux操作与测序原理的实战考核，华大基因2023年9月生物信息分析水平初测试题正是其中的典型代表，整理为面向生信工程师、生物信息分析师求职者的笔试复习资料。测试以单选题、填充题、简答题为主，覆盖两大核心考点：一是Linux文件权限管理，重点考查chmod命令的数字与符号表示法、常见权限状态变化；二是基因组学基础，围绕Sanger测序的双脱氧链终止原理、高准确度优势、成本与速度局限，以及其在验证测序结果、特定基因变异分析中的典型应用。资源压缩包共1个PDF文件，大小3.52MB，部分题目附带个人解析，便于对照考点快速查漏补缺。已有266人学习浏览，适合在面试前系统梳理基础知识点，尤其有助于短期内强化Linux权限操作与经典测序技术的记忆。

1. 生物信息分析水平初测到底在考什么：不是算法，是流程意识

2023年9月那套生物信息分析水平初测试题，我见过很多人倒在第一轮：不是算法不会，而是连一个干净的BAM文件都给不出来。这套题说到底是让你在三到四个小时内，把一份真实的测序数据从FASTQ跑到VCF，再把结果讲清楚。它考的不是你背了多少工具，而是你能不能按工程方式把流程跑完、跑对、跑得可复现。题里不会让你手推动态规划，也不会让你默写BWT的构造过程，但会给你一份可能掺着接头和低质量碱基的真实数据，看你能不能交出一个能用的结果。适合谁？刚转行想做生信的人、实验室里想系统补一遍分析细节的人、以及面试前需要给自己摸底的人。下面这套拆解，按我自己的答题顺序来讲。

2. 初测的四个考核面：从FASTQ质控到变异检测的选型逻辑

2.1 为什么初测永远绕不开FASTQ质控：质量值、接头与FastQC的判读

FASTQ每一行记录一条read，四行一组：序列标识、碱基序列、分隔符、质量字符串。质量字符串每个字符对应一个Phred质量值，ASCII码减去33得到Q值。Q30意味着碱基错误率约千分之一，Q20是百分之一。初测给的数据通常不会故意做成低质量，但真实数据里必然有接头残留、末端质量衰减和重复序列，这三样东西不处理干净，后面的比对和变异检测都会被污染。

最常见的做法是先跑FastQC看一眼全局质量分布。拿到题第一件事不是写脚本，而是先做质控报告。命令很简单：

fastqc sample_R1.fastq.gz sample_R2.fastq.gz -o fastqc_out -t 4

这条命令会为每个输入文件生成一个HTML报告和zip包。-o指定输出目录，-t是并行线程数。跑完先看Per base sequence quality的箱线图，再看Adapter Content和Overrepresented sequences两个模块。如果Adapter Content在read末端明显翘起来，说明接头去除不彻底；如果Overrepresented sequences里有大量相同序列，可能是PCR重复或接头二聚体。

判断完质控报告后，下一步就是清洗。清洗用Trimmomatic比较常见，双端测序的写法：

trimmomatic PE -phred33 \ sample_R1.fastq.gz sample_R2.fastq.gz \ clean_R1.fastq.gz clean_R1_unpaired.fastq.gz \ clean_R2.fastq.gz clean_R2_unpaired.fastq.gz \ ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \ LEADING:3 TRAILING:3 \ SLIDINGWINDOW:4:15 MINLEN:36

-phred33指定质量编码方式，现在绝大多数测序平台都是Phred33，老数据的Phred64偶尔还会在旧题里出现，跑FastQC时如果质量值分布图整体异常偏低，先怀疑这个。ILLUMINACLIP后面的三个数字依次是种子错配数、palindrome模式阈值和simple模式阈值，一般用2:30:10就行。SLIDINGWINDOW:4:15表示从read开头以4碱基为窗口滑动，窗口平均质量低于Q15就从窗口处截断。MINLEN:36是保留最小长度，短于36的read直接丢弃。清洗完再跑一次FastQC，确认Adapter Content已经被压到基线附近再进入比对。

这套质控流程的价值在于：初测的评分点往往不在你最后VCF里有多少个变异，而在你的中间文件是否合理。如果比对率低、插入片段分布异常，大概率是质控环节出了问题。所以拿到数据先质控，不要急着跑比对。

2.2 比对工具怎么选：BWA-MEM与Bowtie2的边界

DNA测序数据的比对，行业内默认首选是BWA-MEM。原因很直接：它容忍错配和indel的能力强，对长读长和短读长都适用，而且速度和内存占用在常规服务器上表现稳定。Bowtie2不是不能用，它的局部比对模式在RNA-seq里有优势，但对DNA重测序来说，BWA-MEM处理indel的敏感度更高，在变异检测场景里更可靠。

先建索引，再比对。很多人在这一步翻车，是因为拿到参考基因组就直接跑bwa mem，忘记先建索引，或者索引和参考基因组版本不匹配。正确顺序：

bwa index ref.fa bwa mem -t 8 \ -R '@RG\tID:sample1\tSM:sample1\tPL:ILLUMINA' \ ref.fa clean_R1.fastq.gz clean_R2.fastq.gz \ | samtools sort -@ 4 -o sample1.sorted.bam -

-R参数是所有比对命令里最不该省的一项，它给read group打标签，后面的GATK流程靠它区分样本。ID是read group唯一标识，SM是样本名，PL是测序平台。多个样本合并分析时，如果没有RG标签，GATK会直接报错或把所有read当成一个样本处理。管道右边直接接samtools sort，省一步生成未排序BAM的中间文件，磁盘占用也能少一些。注意最后的连字符表示从标准输入读取。

samtools sort默认按坐标排序，排序完的BAM才能被GATK识别。比对完成后先做一次快速验证：

samtools flagstat sample1.sorted.bam

看总比对率和properly paired的比例。正常人类WES数据比对率应该在99%以上，properly paired在95%以上。如果比对率掉到90%以下，别急着往下走，先查参考基因组是不是选错了、接头是不是没去干净、或者样本本身有污染。

2.3 变异检测的必考写法：GATK流程与参数为什么要这样设

初测的变异检测环节，几乎绕不开GATK的HaplotypeCaller。它不是唯一的工具，但它是考核里最常出现的标准答案。HaplotypeCaller和旧式UG（UnifiedGenotyper）最大的区别是：它先在每个样本上做局部de novo组装，把活跃区域的read重组成单倍型，再比对参考基因组确定变异位点。这个策略对indel的检测能力明显更强，尤其是在真实数据存在碱基错误和比对模糊时。

GATK4之后的推荐做法是先输出GVCF，再联合批量样本做GenotypeGVCFs。即使初测只给一个样本，也建议按GVCF模式跑：

gatk HaplotypeCaller \ -R ref.fa \ -I sample1.sorted.bam \ -O sample1.g.vcf.gz \ -ERC GVCF \ --native-pair-hmm-threads 8

-ERC GVCF表示输出每个位点的信息，包括所有参考位点和变异位点，后面用GenotypeGVCFs联合时直接使用。如果只跑单样本，也先输出GVCF再转VCF，这样做的好处是中间文件可以复用，后面追加样本时不需要重跑比对。--native-pair-hmm-threads是控制组装时并行线程数的，按机器核数调，别设成超过物理核数。

拿到原始VCF后，需要做过滤。GATK Best Practices推荐的过滤标准这些年在版本间有调整，但核心指标没变过：QD、FS、MQ、MQRankSum。保守但稳妥的过滤参数：

gatk VariantFiltration \ -R ref.fa \ -V sample1.g.vcf.gz \ -O sample1.filtered.vcf.gz \ --filter-expression "QD < 2.0 || FS > 60.0 || MQ < 40.0" \ --filter-name "LowQual"

过滤表达式里QD是质量值除以深度，indel位点往往偏低；FS是链特异性检验，存在链偏向的变异是假阳性的概率高；MQ是比对质量，低于40说明这个位点的read比对质量整体不行。初测不会要求你把阈值调到多完美，但如果你能说出为什么用这个阈值，比闷头跑完一个流程得分高。过滤只是给位点打上LowQual标签，不是删除位点，后续下游分析时只需要保留PASS位点即可。

3. 把一套典型初测题跑通：从原始数据到VCF的完整命令流

3.1 准备参考基因组与索引：跑题之前先排掉的三个暗坑

初测通常会提前告知参考基因组的物种和版本，但不会帮你准备好索引。下载参考基因组后，第一件事是检查它的染色体命名，这一步能省下后面好几个小时的排查时间。人类参考基因组有两个常见的命名体系：UCSC风格是chr1、chr2、chrX，Ensembl风格是1、2、X。如果参考基因组用的是Ensembl编号，而你后续的bed文件、样本信息表里写的是chr1，所有注释步骤都会对不上。

第二个暗坑是参考基因组文件里可能带decoy序列或者未定位的contig（比如chrUn_gl000220），这些序列不影响主流程，但会让比对和GATK的interval处理变慢。如果初测给了目标区域bed文件，直接用bed把分析限定在目标区间内，不要全基因组跑。

第三个暗坑是bwa index的构建时间。人类参考基因组约3.1Gb，bwa index在普通服务器上可能跑一两个小时。合理的做法是拿到题后第一时间启动索引构建，让它在后台跑着，同时去处理质控和清洗。索引命令：

bwa index -a bwtsw ref.fa samtools faidx ref.fa

-a bwtsw适用于大于2Gb的大基因组，人类参考基因组必须加这个参数，小基因组不加也能跑，但大基因组不加会报错。samtools faidx生成.fai索引文件，后面GATK和samtools view按区间提取时需要它。两个索引都建好后，用下面的命令快速确认：

samtools view -H sample1.sorted.bam | grep "^@SQ" | head -5

看输出里的SN字段是否和参考基因组fasta的染色体名一致。如果在比对完成之后才发现染色体名对不上，那后悔药只有重跑比对，所以这一步检查务必放在比对之前。

3.2 bwa + samtools 的最小比对流程：命令顺序为什么不能乱

比对的完整流程是：建索引、比对、排序、标记重复、建立BAM索引。排序和标记重复的顺序不能颠倒，因为MarkDuplicates需要按坐标排序后的BAM才能正确识别彼此相邻的重复read。只做初测的话，标记重复这步未必强制，很多流程为了省时间会跳过，但如果你交出的BAM里没标记重复，而下一轮分析要用它做变异检测，PCR重复会影响等位基因频率的准确性。

推荐的最简流程是只做排序和建索引，重复标记放在有时间余量时补做：

bwa mem -t 8 -M \ -R '@RG\tID:sample1\tSM:sample1\tPL:ILLUMINA' \ ref.fa clean_R1.fastq.gz clean_R2.fastq.gz \ > aligned.sam samtools sort -@ 4 -o sample1.sorted.bam aligned.sam samtools index sample1.sorted.bam

-M参数会把长的split read标记为secondary alignment，这是GATK后续处理需要的标记方式。先输出SAM再排序和直接管道输出排序结果的区别在于：管道写法省磁盘空间，但排错时不方便看SAM内容。我一般先落盘SAM，确认比对没问题后再排序，只是个人习惯，磁盘紧张时直接管道更稳。

排序后立刻做三次检查，任何一次异常都别往下走。第一看flagstat的properly paired比例，第二看覆盖深度，第三看插入片段分布：

samtools flagstat sample1.sorted.bam samtools depth -r chr1:1000000-1001000 sample1.sorted.bam | awk '{sum+=$3} END {print sum/NR}'

第二条命令只看一个1kb区间，用来快速估算深度。如果这个值异常低或者为0，最大可能是参考基因组染色体命名不一致，或者reads没有比对上去。插入片段分布的检查在QC的FastQC报告里能间接看到，BAM层面的检查用samtools stats：

samtools stats sample1.sorted.bam | grep "insert size average"

插入片段平均长度在300-500bp范围内是正常双端测序的表现。如果平均值只有几十bp，说明文库构建有问题，或者数据本身是片段化的，这个信息会在后续变异检测时影响对indel的解读。

3.3 用GATK HaplotypeCaller出VCF：RG标签与过滤参数一次说清

比对完成后进入变异检测。第一步先对BAM做BaseRecalibrator（BQSR），这一步在纯初测里不是必须，因为BQSR需要已知的高置信度变异集（比如某个物种的dbSNP和gold standard indel），如果题里没提供，强行跑反而会把参考集的不准带进结果。没有已知变异集时，直接跑HaplotypeCaller是更安全的做法。

HaplotypeCaller需要BAM文件头里有正确的read group信息。检查方法：

samtools view -H sample1.sorted.bam | grep "@RG"

如果输出为空或者没有SM标签，需要用samtools addreplacerg补上：

samtools addreplacerg \ -r 'ID:sample1' -r 'SM:sample1' -r 'PL:ILLUMINA' \ -o sample1.rg.bam sample1.sorted.bam samtools index sample1.rg.bam

addreplacerg会生成一个新BAM文件，原文件保留，磁盘空间紧张的话可以加--no-PG参数省掉注释。这道工序做完后重新跑一次samtools view -H确认，再进入GATK。很多人在这一步翻车，GATK运行时报"read group is missing SM tag"错误，浪费半小时。

HaplotypeCaller跑完后，用bcftools做一次快速统计：

bcftools stats sample1.g.vcf.gz > vcf_stats.txt grep "number of SNPs:" vcf_stats.txt grep "number of indels:" vcf_stats.txt

预期SNP数量和indel数量取决于物种和你分析的是全基因组还是全外显子。全外显子人类数据SNP一般在一两万个量级，如果结果只有几十个，先回去查比对率；如果结果有几十万个，先怀疑参考基因组版本是不是和样本物种不匹配。变异数量异常通常不是HaplotypeCaller的问题，而是前面所有步骤累积的偏差。

3.4 用Python做结果统计与可视化：初测里最容易被忽略的得分点

初测通常不只让跑出VCF，还会要求做一点统计描述，比如突变类型分布、转换/颠换比（Ti/Tv）、突变频率分布。这部分用Python可以快速搞定。读取VCF最稳的方式是用pysam，它对VCF/VCF.gz的解析比纯文本逐行读要可靠得多：

import pysam from collections import Counter vcffile = pysam.VariantFile("sample1.filtered.vcf.gz") ti_tv_counter = Counter() type_counter = Counter() for rec in vcffile.fetch(): if rec.filter.keys() != ["PASS"]: continue ref = rec.ref alt = rec.alts[0] if len(ref) == 1 and len(alt) == 1: type_counter["SNP"] += 1 transition = {("A", "G"), ("G", "A"), ("C", "T"), ("T", "C")} if (ref, alt) in transition: ti_tv_counter["Ti"] += 1 else: ti_tv_counter["Tv"] += 1 else: type_counter["Indel"] += 1 print("SNP:", type_counter["SNP"], "Indel:", type_counter["Indel"]) print("Ti/Tv:", ti_tv_counter["Ti"] / max(ti_tv_counter["Tv"], 1))

脚本先过滤掉非PASS位点，再按ref和alt的长度区分SNP与indel，最后统计转换和颠换。人类全外显子数据的Ti/Tv比值通常在2.0到3.0之间，低于2.0说明可能有假阳性混杂，高于3.5说明过滤过严或数据有偏差。这个数字在初测评分里经常被当作结果质量的一票否决项。

如果题里要求画图，用matplotlib画一个变异密度分布，按染色体位置分bin统计：

import matplotlib.pyplot as plt bins = {} for rec in vcffile.fetch(): if rec.filter.keys() != ["PASS"]: continue chrom = rec.chrom pos = rec.pos if chrom not in bins: bins[chrom] = [] bins[chrom].append(pos) fig, ax = plt.subplots(figsize=(10, 6)) for chrom, positions in bins.items(): ax.hist(positions, bins=100, alpha=0.6, label=str(chrom)) ax.set_xlabel("Genomic Position") ax.set_ylabel("Variant Count") ax.legend() plt.tight_layout() plt.savefig("variant_density.png", dpi=150)

这个图能看到是否有区域异常聚集的变异hotspot，初测数据如果来自真实样本，X染色体和线粒体的变异密度会和常染色体不同，这属于正常生物学信号，不用慌。

4. 初测避坑指南：五条样本级别的血泪经验

4.1 FastQC接头污染爆红，但实际是剪不干净

现象：FastQC的Adapter Content模块在read末端明显翘起，Overrepresented sequences里有大量接头序列。清洗后再跑FastQC，接头比例依然超标。

原因：Trimmomatic的ILLUMINACLIP参数设置得太宽松，或者接头文件里没有包含正确的接头序列。TruSeq3-PE.fa只覆盖了Illumina标准接头，如果建库试剂盒用的是定制接头，这个文件匹配不上。另外palindrome模式阈值设太高也会导致接头未被识别。

解决：先看FastQC报告里Overrepresented sequences的具体序列是什么，比对到TruSeq3-PE.fa里的序列，确认是否匹配。如果序列不一致，用单端模式手动截掉接头序列，或者调整ILLUMINACLIP参数里的种子错配数从2降到1。更快的办法是用cutadapt重新清洗，直接把已知接头序列传给它的-a和-b参数。

4.2 比对率低到怀疑数据：参考基因组版本与染色体命名

现象：samtools flagstat显示比对率只有70%-85%，properly paired比例更低，但FastQC的碱基质量看着正常。

原因：参考基因组下载错了版本，或者用的是同一个物种但不同品系的基因组。人类参考基因组从GRCh37到GRCh38，虽然整体差异不大，但某些区域的序列结构完全不同，比对时会成片map不上。还有一种情况是参考基因组用的是ENSEMBL编号，而样本的bed文件或注释文件里全是chr前缀，导致GATK的interval处理时找不到对应染色体。

解决：确认参考基因组的准确来源，重新下载与题面对应版本的fasta。用samtools view -H看BAM的@SQ行，和参考基因组的序列名逐一对比。如果只是命名体系不一致，可以用samtools reheader重写头部替换染色体名，但强烈不建议这么做，至少我不推荐——重写头后序列顺序和坐标可能产生偏移，还不如重新比对来的干净。

4.3 HaplotypeCaller空VCF：RG标签缺失的经典翻车

现象：HaplotypeCaller跑完没有任何报错，但输出的VCF里只有几行记录，甚至一个变异都没有。

原因：BAM文件头里缺少read group信息，GATK在底层处理时把不同的read当成了多个样本的混合数据，内部会在分组后丢弃几乎所有位点。这个报错在GATK3时代会明确提示，GATK4有时只给warning，导致流程继续跑完却得到空结果。

解决：比对之前就给BAM加上RG标签，不要等到HaplotypeCaller报错才想起来。已经跑完的话，用samtools addreplacerg补上标签后重新建索引，再跑一次HaplotypeCaller。补标签后一定要重新生成bai索引，否则GATK读取时还是会报错。这个坑我在真实的某次项目里踩过一遍，空VCF加排查日志，前前后后浪费了一个下午。

4.4 Python读VCF内存爆炸：流式处理与格式先行

现象：用pandas.read_csv直接读VCF文件，几十G的VCF.gz读进内存后进程被杀，或者耗时极长。初测数据量不大时可能不触发，但遇到全基因组级别的VCF就必然出事。

原因：VCF是行列式文本，全基因组VCF有几千万行，pandas会把所有行读进内存，结构变异注释列又多，内存直接被打满。用纯文本逐行解析也有同样问题，因为数据量太大时逐行读虽然省内存，但速度慢到无法接受。

解决：用pysam或cyvcf2这种专为VCF设计的解析库，流式迭代处理，不一次性加载全文件。上面3.4节的写法就是流式的，每处理完一行就释放一行，不会累积。如果只是要统计位点数量，用bcftools stats更快，不需要写Python。记住一条原则：能用bcftools完成的统计不要写脚本，写了脚本也要先检查有没有现成工具可用。

4.5 结果被判定分析失败：可复现性也是评分项

现象：结果文件都齐了，VCF数量也正常，但评审反馈说分析不可靠，要求重做。

原因：没有记录用的什么版本的工具、什么参数、参考基因组的下载地址和md5值。初测看的不只是最终结果，还看过程的可复现性。如果下次换台机器重新跑，连你自己都说不清当时用的哪个版本的GATK，那结果就是黑匣子里出来的，不具备可信度。

解决：从开工第一秒就养成分步记录的习惯。我的做法是在工作目录下放一个commands.sh文件，每执行一步就把命令追加进去，并在关键步骤后面写注释说明版本和参数含义。所有中间文件保留，不要随手删除。最后交付时附上一个简短的README，写清楚每个文件是哪个环节产生的。这个习惯在做大项目时也能救命，某次某公司一个真实项目因为没记录参数，复现的时候发现GATK版本差异导致结果对不上，折腾了三天才排查清楚。

5. 把初测当工程做：用日志和中间文件反向验证结果

5.1 用 set -euo pipefail 和日志文件守住流程底线

写流程脚本时，第一行就加上set -euo pipefail。这条命令的意思是：任何一条命令返回非零状态码，脚本立即退出；使用未定义的变量直接报错；管道中任何一条命令失败，整个管道的退出码就是失败的。没有这三条，脚本会在比对失败后继续往下跑，最后用残缺的BAM文件生成VCF，结果全错还找不到原因。配合日志输出的做法是给每个步骤单独写日志：

set -euo pipefail exec 1> >(tee run.log) 2>&1

tee在这里把标准输出同时打到屏幕和日志文件，排查问题时直接看run.log的时间戳和命令顺序。每跑完一个步骤，在日志里打印一行标记，比如"QC done, 3 files passed"，这样失败时可以快速定位到具体是哪一步。

5.2 三个必看的中间验证点：flagstat、stats与随机抽样

我习惯在三个节点停下验证，不闷头跑完全流程。比对完马上看flagstat的properly paired比例，低于95%就回头查参考基因组和清洗参数。HaplotypeCaller跑完看bcftools stats的Ti/Tv比值，低于2.0就检查过滤参数是不是过松。最后在交付前随机抽一个区域，用IGV可视化看一下真实read堆积情况，确认变异位点确实有read支持。

samtools view -b sample1.sorted.bam chr1:10000000-10001000 | samtools depth - | awk '{sum+=$3} END {print "depth:", sum/NR}'

这条命令随机抽一个10kb区域算深度，如果深度为0或异常低，说明这个区域在比对时被整体丢弃了，需要检查是不是目标区域bed文件里的坐标和参考基因组命名不匹配。这步检查做完，结果才有底气交付。

以前我在做某次模拟项目X时，就是因为没加set -euo pipefail，某个样本比对失败后脚本继续往后跑，最后VCF全是空的，白等一整夜。后来养成了每个脚本开头先写这一条，现在闭着眼都能打出来。这个习惯帮我少走了很多弯路，也希望帮到你。

本文还有配套的精品资源，点击获取