序列比对与处理：BWA-MEM比对参数选择，SAM/BAM格式，排序，标记重复
序列比对这一步，说白了就是要把你手里那些乱七八糟的测序片段，老老实实贴到参考基因组上去。很多新手以为只要跑个bwa mem就完事了，结果后面的变异检测全是坑。我刚开始接触这个的时候，也被参数调教得够呛。
1. BWA-MEM 核心参数详解
BWA-MEM是目前最常用的比对工具，专门针对70bp以上的测序数据。我个人习惯先跑一遍默认参数看看效果，但真正上线跑项目的时候，一定会针对数据特征做微调。
| 参数 | 默认值 | 作用场景 |
|---|---|---|
| -k | 19 | 种子长度，短读长（<100bp）建议调低至17-18 |
| -r | 1.5 | 链结惩罚分，重复区域较高的数据可适当提高 |
| -c | 10000 | 对高覆盖度数据可降至5000-8000，减少假阳性 |
| -M | 关闭 | 将split read标记为副比对，做SV检测时建议开启 |
| -t | 1 | 线程数，根据CPU核心数设定，建议不超过24 |
大家在看BWA参数时，别光盯着默认值。我记得有次客户给的WGS数据，初始比对率只有94%，排查后发现就是-k这个种子长度问题。样本测序读长只有76bp，种子长度默认19，导致大量短片段无法锚定。调成17后，比对率直接拉到98%。
常用比对命令模板
bwa mem -t 16 -k 17 -M -R '@RG\tID:样本A\tSM:样本A\tLB:文库1\tPL:ILLUMINA' \
ref.fasta sample_R1.fastq.gz sample_R2.fastq.gz > aligned.sam这里面最容易被忽略的就是RG标签。嗯，这里要注意，RG标签的SM字段决定了后续样本名称，如果多个样本用了同一个SM值，下游工具会认为它们是同一个样本。我曾经因为这个问题把一个队列的结果全部搞混，重跑了一整天。
2. SAM与BAM格式基础
比对完出来的SAM文件，说白了就是一个巨大的文本表格，每一行代表一条比对记录。它的结构其实很清晰，但第一次打开时可能会被吓到。
SAM格式的每一行主要包含11个必填字段：
- QNAME - 测序片段ID，来源fastq的read name
- FLAG - 位标记，告诉你比对状态
- RNAME - 参考序列ID
- POS - 比对到参考的起始位置
- MAPQ - 比对质量值
- CIGAR - 比对操作符码
- RNEXT - mate pair的参考序列
- PNEXT - mate pair的起始位置
- TLEN - 模板长度
- SEQ - 原始序列
- QUAL - 测序质量值字符串
FLAG字段是很多初学者最头疼的地方。实际上它是一个二进制位图，用总和来表示比对状态。举个简单的例子，FLAG=99表示这条read是paired-end的第一次读段，并且与参考序列正链比对。你想想看，99 = 1（paired）+ 2（proper pair）+ 32（未比对）+ 64（第一次读段）= 99，就是这么算的。
你可以在终端里直接用 samtools flags 命令查看FLAG含义。比如输入 samtools flags 99，它会告诉你这条记录的完整状态。
BAM是SAM的二进制压缩格式，占用空间只有SAM的1/5左右。几乎所有下游分析工具都要输入BAM格式。我自己的流程中从来不会保留SAM文件，比对一结束就转成BAM并删除原始SAM。
3. Samtools排序：为什么一定要做？
很多人觉得排序是个麻烦事，但这事非做不可。BWA输出的比对结果默认是按读段ID排列的，而后续的变异检测工具（比如GATK、FreeBayes）需要按基因组坐标排序才能正常工作。为什么？你想想看，按坐标排序后，相邻的比对记录在物理位置上也是相邻的，这样才能高效地扫描基因组上的每一个位点。
排序的命令很简单：
samtools sort -@ 8 -m 2G -o sorted.bam aligned.sam参数说明：
-@ 8使用8个线程，加快排序速度-m 2G每个线程分配2GB内存，总和不要超过物理内存-o输出文件路径
我曾尝试过对高深度WGS数据（约60X）用单线程排序，结果跑了两个小时还没完。调整成8线程、每线程分配4G内存，十来分钟就搞定了。内存和线程数是排序的瓶颈，一定要合理分配。
注意：如果内存不足，samtools sort会创建大量临时文件，严重拖慢速度。建议临时目录（通过 -T 参数指定）放在SSD上，普通机械盘会慢到让你怀疑人生。
4. 标记重复（MarkDuplicates）
标记重复是变异检测流程中非常关键的一步，但很多人不理解它的意义。
PCR重复，说白了就是建库过程中同一个DNA片段被多次扩增，导致测序数据中出现多条完全相同的read。这些重复read会显著影响变异检测的准确性，尤其是会高估低频变异频率。
Picard工具集的MarkDuplicates是目前最常用的标记重复工具：
picard MarkDuplicates \
I=sorted.bam \
O=marked_dup.bam \
M=dup_metrics.txt \
REMOVE_DUPLICATES=false \
ASSUME_SORT_ORDER=coordinate注意上面的 REMOVE_DUPLICATES=false，
为什么我建议不要直接删除重复？
- 保留重复但标记，后续变异检测工具可以自己决定是否忽略这些位点
- 直接删除后如需溯源，原始数据已经丢失
- 某些特殊应用（如CNV分析）需要保留重复信息
我曾经接手过一个项目，之前的生信人员直接设置REMOVE_DUPLICATES=true，把PCR重复全删了。后来做低频变异验证时发现很多位点深度不够，再想恢复原始比对文件已经来不及了。从那以后我的流程里永远只标记、不删除。
标记重复的原理并不复杂。工具会比对read pairs的5'端坐标，如果两个read pair的5'坐标完全一致，就认为是PCR重复。真正的生物学重复（即同一位置上独立的两个DNA片段）概率极低，可以忽略不计。
本章知识图谱
下面这张图展示了从原始测序数据到变异检测之间的处理流程：
这张图的路径也很简单：原始数据先经过BWA-MEM比对生成SAM，转成BAM格式后必须先做排序，然后才能标记重复。有朋友可能想问，能不能跳过排序直接标记重复？答案是不行。MarkDuplicates要求输入数据按坐标排序过，否则它无法判断两条read是否是重复。
常见参数调优经验
针对不同数据类型的建议
| 数据类型 | BWA参数调整 | 排序策略 | 重复标记策略 |
|---|---|---|---|
| WGS (30X) | -k 19, -r 1.5, -M | 8线程, 2G/线程 | 保留标记，不删除 |
| WES (100X+) | -k 17, -r 2.0, -M | 4线程, 1G/线程 | 保留标记，但重复率过高时考虑删除 |
| RNA-seq | 建议使用STAR/HISAT2 | 8线程, 2G/线程 | 通常不标记重复 |
| 超深度靶向测序 (1000X+) | -k 15, -c 5000 | 16线程, 4G/线程 | 必须标记，建议删除重复 |
看到这里你应该能感受到，比对和处理流程虽然看起来简单，但每步都有讲究。我的建议是，先用小数据集（比如一个样本的1M条read）跑一遍完整流程，用samtools stats检查比对率和重复率，确认参数设置没问题再放量跑全样本。血的教训告诉我，批量跑完才发现参数有问题的感觉，谁都不想来第二回。