十年匠心定制 · 商业建站与技术教学双线并行 咨询热线:400-886-1026 service@lmnt.cn
ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

Fastp质控实战:从参数调优到批量处理的完整指南

Fastp质控实战:从参数调优到批量处理的完整指南 测序数据下机之后第一件事永远是看质量。不管你是做全基因组、转录组还是扩增子拿到手的 fastq 文件里藏着接头污染、低质量碱基、长度分布异常这些麻烦直接拿去比对或者组装结果大概率会让你怀疑人生。Fastp 就是在这个环节里被大量使用的工具之一它把质控、过滤、接头去除、UMI 处理这些步骤揉进了一个命令里速度还快得离谱。这篇内容面向的是刚接触 fastq 质控、或者已经在用 Fastp 但想把它用透的人我会从实际项目里的操作链路出发把参数选择、UMI 处理、批量脚本、结果解读这些环节拆开讲清楚尽量让你看完就能直接抄作业。1. 先搞清楚 Fastp 到底替你干了哪些活很多人对 Fastp 的印象停留在一个质控工具但真到项目里你会发现它承担的角色比想象中多。理解它做了什么才能知道哪些参数值得调、哪些默认值可以直接用。1.1 从 fastq 到 clean reads 的完整链路一个典型的 fastq 文件里每条 read 由四行组成序列标识行、碱基序列、正号分隔行、质量值行。质量值用 ASCII 字符表示常见的是 Phred33 编码字符!对应质量值 0I对应 40。Fastp 读入这些数据后会依次做几件事先扫描全局质量分布再逐条 read 做滑窗质量修剪然后识别并切除接头序列接着按长度和复杂度过滤最后输出 clean fastq 和一份 HTML 报告。这个顺序不是随便定的。滑窗修剪放在接头切除之前是因为接头区域的碱基质量往往偏低先修剪能减少接头识别时的干扰。而长度过滤放在最后是因为前面的修剪和切除都会改变 read 长度如果提前过滤可能把本来能救回来的 read 误杀。理解这个执行顺序你在调参时就不会把参数放错位置。Fastp 的滑窗逻辑默认是 4 个碱基一个窗口从 read 的 5 端往 3 端滑动计算窗口内平均质量低于阈值就把窗口及之后的部分切掉。这个设计对 Illumina 数据特别友好因为它的质量下降通常集中在 3 端。我实测过一批双端 150bp 的数据开启滑窗后 Q30 比例从 88% 提到了 94%代价是平均读长掉了大约 12bp这个取舍在大多数分析里是划算的。1.2 为什么它比传统组合拳更省事以前做质控常见做法是 Trimmomatic 去接头加 FastQC 看质量再写脚本统计过滤前后数据量。这套流程能跑但有两个痛点一是中间文件多磁盘占用大二是 FastQC 只给报告不给过滤能力你得自己判断阈值再回头调 Trimmomatic 参数来回折腾。Fastp 把这两件事合并了。它一边过滤一边生成 HTML 报告报告里直接给出过滤前后的质量曲线、碱基组成、长度分布、接头统计。你跑完一条命令既拿到了 clean data又拿到了判断依据。更关键的是它支持多线程-w参数指定线程数后处理 30GB 的双端数据在 16 核机器上大概十几分钟就能跑完这个速度在批量项目里能省下大量等待时间。还有一个容易被忽略的点Fastp 默认会自动检测接头序列不需要你手动提供 adapter 列表。它通过比对 read 之间的重叠区域来推断接头对标准 Illumina 文库的通用接头识别率很高。当然如果你的文库类型特殊比如用了定制接头那就得用--adapter_sequence手动指定这个后面会细说。1.3 适用场景与不适用场景Fastp 最适合的是 Illumina 短读长数据的常规质控包括 WGS、WES、RNA-seq、ChIP-seq、扩增子这些。它对双端数据的处理尤其成熟能利用 read overlap 做接头检测和纠错。但它不是万能的。长读长数据比如 PacBio 和 NanoporeFastp 的处理逻辑就不太匹配这类数据更适合用专门的工具。另外如果你的需求是做严格的去重或者复杂的 UMI 纠错Fastp 的 UMI 功能是够用的但如果是高度复杂的 UMI 设计可能需要配合其他工具一起用。还有一个边界Fastp 不做比对所以它无法识别污染物种那是比对工具的活。2. 参数不是越多越好关键参数逐个拆Fastp 的参数列表很长但真正常用的就那么十几个。我把它们分成质量过滤、接头处理、长度过滤、UMI 处理四组来讲每个参数说清楚调它会发生什么。2.1 质量过滤参数滑窗与碱基阈值-q或--qualified_quality_phred控制的是碱基质量阈值默认 15。意思是质量值低于 15 的碱基会被标记为不合格。-u或--unqualified_percent_limit默认 40表示一条 read 里不合格碱基超过 40% 就整条丢弃。这两个参数配合使用。如果你把-q提到 20不合格碱基的定义变严同样的-u下会有更多 read 被丢弃。我一般建议-q保持 15 到 20 之间-u保持 40 左右除非你的下游分析对错误率极其敏感。滑窗相关的参数是-W窗口大小默认 4和-M窗口平均质量阈值默认 20。-M是滑窗修剪的核心它决定从哪个位置开始切。如果你发现 clean data 的 3 端质量还是不理想可以把-M提到 25 甚至 28但要注意读长会进一步缩短。我做过一组对比-M 20时平均读长 138bp-M 25时降到 131bpQ30 从 93% 升到 96%。这个取舍取决于你的下游是比对还是组装比对通常能容忍短一点但质量更高的 read。还有一个参数-e或--average_qual默认 0 表示不启用。它按整条 read 的平均质量过滤如果你设成 20平均质量低于 20 的 read 直接丢。这个参数适合在数据质量整体偏差时做粗筛但设太高会损失大量数据慎用。2.2 接头处理自动检测与手动指定Fastp 默认开启接头检测靠的是--detect_adapter_for_pe这个行为双端数据自动启用。它的原理是找 read1 和 read2 之间的重叠区域如果重叠区之外还有序列那多半就是接头。这个方法对标准文库很有效但对插入片段很短的文库可能误判。手动指定接头用--adapter_sequenceread1 接头和--adapter_sequence_r2read2 接头。Illumina 通用接头序列是AGATCGGAAGAGCACACGTCTGAACTCCAGTCAread1和AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTread2。如果你不确定文库用的什么接头可以先跑一次默认检测看 HTML 报告里的接头统计如果检测到的接头比例异常高或者序列奇怪再考虑手动指定。接头切除还有几个辅助参数。--adapter_fasta可以提供一个 fasta 文件里面放多个接头序列适合多重接头的场景。-x或--trim_poly_x用于去除 polyX 尾比如 polyA 尾这个在 RNA-seq 里比较常见。--cut_front、--cut_tail、--cut_right是三个滑动切割参数分别从 5 端、3 端、或者向右滑动切割一般用默认关闭就行除非你有明确的切割位置需求。2.3 长度与复杂度过滤别把好数据误杀了-l或--length_required默认 15表示修剪后长度低于 15 的 read 会被丢弃。这个值设太低会留下很多碎片设太高会损失数据。对于 150bp 双端数据我一般设 30 到 36因为低于 30bp 的 read 在比对时基本无法唯一定位。--length_limit默认 0 表示不限制如果你设了值超过这个长度的 read 会被丢弃。这个参数用得少主要针对异常长的 read。复杂度过滤用-Y或--low_complexity_filter默认关闭。开启后复杂度低于-y默认 30%的 read 会被过滤。复杂度计算方式是相邻碱基不同的比例比如AAAAAA的复杂度是 0ATATAT的复杂度是 100%。低复杂度序列在比对时容易产生多重映射如果你做的是 SNP calling 或者表达定量开启这个过滤能减少假阳性。但要注意某些真实序列本身复杂度就低比如 polyA 尾或者简单的重复序列开启后可能误伤。2.4 UMI 处理从理论到参数落地UMI 是 Unique Molecular Identifier一段随机序列用来标记原始 DNA 分子主要解决 PCR 扩增偏倚和去重问题。Fastp 支持 UMI 处理但前提是你的 UMI 位置和结构得先搞清楚。常见的 UMI 设计有两种一种是 UMI 在 read 的 5 端比如UMI 插入片段另一种是 UMI 在 index 里通过--umi_loc指定位置。Fastp 的 UMI 参数包括--umi启用 UMI 处理、--umi_locUMI 位置可选 index1、index2、read1、read2、per_index、per_read、--umi_lenUMI 长度、--umi_prefixUMI 前缀用于输出标识。举个例子如果你的 UMI 是 read1 的前 8bp命令里加--umi --umi_locread1 --umi_len8。Fastp 会把这段 UMI 切下来移到 read 名称里格式类似UMI_8bp_序列。这样后续去重工具就能根据 UMI 识别同一分子的不同 read。这里有个坑UMI 处理必须在接头切除之前完成因为 UMI 就在 read 末端如果先切接头可能把 UMI 一起切掉。Fastp 内部会处理这个顺序但你要确保--umi参数正确否则 UMI 信息会丢失。另外如果 UMI 在 index 里你需要先用--umi_locindex1或index2并且确保 index 读取正确。3. 单端与双端、单样本与批量的命令组织参数懂了接下来是怎么把它们组织成能跑的命令。单样本和批量、单端和双端命令结构有差异我分开说。3.1 双端数据的基础命令模板双端是最常见的情况基础命令长这样fastp -i sample_R1.fastq.gz -I sample_R2.fastq.gz \ -o sample_R1.clean.fastq.gz -O sample_R2.clean.fastq.gz \ -w 8 \ -q 20 -u 40 \ -l 36 \ -M 25 \ --detect_adapter_for_pe \ -j sample.fastp.json -h sample.fastp.html逐段解释-i和-I是输入 read1 和 read2-o和-O是输出。-w 8用 8 个线程。-q 20 -u 40是质量过滤。-l 36是长度过滤。-M 25是滑窗质量阈值。--detect_adapter_for_pe显式开启双端接头检测其实双端默认就开写出来更清楚。-j和-h分别输出 JSON 和 HTML 报告。这里有个细节输入输出都用.gz压缩格式Fastp 会自动识别并处理不需要额外解压。输出压缩能省不少磁盘代价是稍微多花一点 CPU 时间但相比磁盘 IO 的节省这个代价可以忽略。3.2 单端数据的参数差异单端数据没有 read2命令简化成fastp -i sample.fastq.gz \ -o sample.clean.fastq.gz \ -w 8 \ -q 20 -u 40 \ -l 36 \ -M 25 \ --adapter_sequence AGATCGGAAGAGCACACGTCTGAACTCCAGTCA \ -j sample.fastp.json -h sample.fastp.html单端没有 read overlap所以自动接头检测不适用必须手动指定--adapter_sequence。如果你不知道接头序列可以先不指定跑一次看报告里的接头统计或者用其他工具先检测。单端的滑窗修剪逻辑和双端一样但因为没有配对信息长度过滤要更谨慎设太低会留下很多无法比对的短 read。3.3 批量处理的 Shell 脚本写法项目里通常有几十上百个样本一个个跑不现实。用 Shell 脚本批量处理是标准做法。下面这个脚本遍历当前目录下所有_R1.fastq.gz文件自动配对 R2输出到clean目录#!/bin/bash mkdir -p clean mkdir -p reports for r1 in *_R1.fastq.gz; do sample${r1%_R1.fastq.gz} r2${sample}_R2.fastq.gz if [ ! -f $r2 ]; then echo 警告$r2 不存在跳过 $sample continue fi echo 正在处理 $sample ... fastp -i $r1 -I $r2 \ -o clean/${sample}_R1.clean.fastq.gz \ -O clean/${sample}_R2.clean.fastq.gz \ -w 8 \ -q 20 -u 40 \ -l 36 \ -M 25 \ --detect_adapter_for_pe \ -j reports/${sample}.fastp.json \ -h reports/${sample}.fastp.html \ 2 reports/${sample}.fastp.log if [ $? -eq 0 ]; then echo $sample 处理完成 else echo $sample 处理失败请检查日志 fi done echo 全部样本处理结束这个脚本里有几个实用点。${r1%_R1.fastq.gz}是 Shell 的参数扩展去掉后缀得到样本名比用sed或basename更高效。2把标准错误重定向到日志文件方便排查问题。$?检查上一条命令的退出状态判断成功与否。如果你样本量特别大串行跑太慢可以用xargs并行ls *_R1.fastq.gz | sed s/_R1.fastq.gz// | \ xargs -I {} -P 4 bash -c fastp -i {}_R1.fastq.gz -I {}_R2.fastq.gz \ -o clean/{}_R1.clean.fastq.gz -O clean/{}_R2.clean.fastq.gz \ -w 4 -q 20 -u 40 -l 36 -M 25 \ --detect_adapter_for_pe \ -j reports/{}.fastp.json -h reports/{}.fastp.html -P 4表示同时跑 4 个样本每个样本用 4 个线程。总线程数是 16和机器核数匹配。这里要注意并行数乘以每样本线程数不要超过机器总核数否则会互相抢资源反而变慢。4. 跑完之后怎么看结果、怎么判断质控是否合格命令跑完不是结束报告里的信息才是决定这批数据能不能用的关键。Fastp 的 HTML 报告信息量很大我挑几个最需要关注的模块讲。4.1 HTML 报告里的关键指标打开 HTML 报告最上面是 Summary 部分直接给出过滤前后的 reads 数、碱基数和 Q30 比例。我一般先看三个数过滤后 reads 保留率、Q30 比例、平均读长。保留率在 80% 到 95% 之间通常算正常。如果低于 70%说明过滤太狠或者原始数据质量差需要回头检查参数。Q30 比例在过滤后应该明显提升如果提升不明显可能是滑窗阈值设太低。平均读长如果掉得太多比如从 150bp 掉到 100bp 以下说明修剪过度要考虑放宽-M。往下看是 Before filtering 和 After filtering 的对比图包括质量曲线、碱基组成、长度分布。质量曲线如果过滤后在 3 端仍然明显下降说明滑窗没切干净可以加大-M。碱基组成图里如果某个碱基在特定位置异常高可能是接头残留或者污染。长度分布图能看出插入片段大小如果分布很窄或者有异常峰可能文库构建有问题。4.2 接头检测结果的解读报告里有一个 Adapter 模块显示检测到的接头序列和比例。如果接头比例超过 5%说明接头污染比较严重需要确认接头切除是否生效。如果检测到的接头序列和你预期的通用接头不一致可能是文库用了定制接头需要手动指定。有个常见情况双端数据里 read1 和 read2 的接头比例差异很大。这通常是因为插入片段长度分布不均短插入片段的 read 更容易读到接头。这种情况下自动检测一般能处理但如果差异极端可以考虑用--adapter_sequence和--adapter_sequence_r2分别指定。4.3 从 JSON 报告里提取批量统计HTML 报告适合单个样本查看批量项目里你需要从 JSON 报告里提取数据做汇总。Fastp 的 JSON 报告结构清晰可以用jq工具解析。下面这个命令提取每个样本的过滤前后 reads 数和 Q30 比例for json in reports/*.fastp.json; do sample$(basename $json .fastp.json) before_reads$(jq .summary.before_filtering.total_reads $json) after_reads$(jq .summary.after_filtering.total_reads $json) q30_before$(jq .summary.before_filtering.q30_rate $json) q30_after$(jq .summary.after_filtering.q30_rate $json) echo -e ${sample}\t${before_reads}\t${after_reads}\t${q30_before}\t${q30_after} done summary.tsv这个汇总表可以直接导入 Excel 或者 R 做可视化。我习惯再加一列计算保留率方便快速筛出异常样本。如果某个样本保留率明显低于其他样本就要单独看它的 HTML 报告找原因。5. 那些文档里不写、但实际会踩的坑参数和命令都清楚了但实际跑起来还是会遇到各种意外。这一节我把自己踩过的坑和排查思路整理出来希望能帮你少走弯路。5.1 内存与线程的隐性冲突Fastp 的内存占用和线程数、数据量都相关。我遇到过在 32GB 内存的机器上跑 8 线程处理 50GB 双端数据时被 OOM killer 杀掉的情况。原因是 Fastp 在接头检测阶段需要缓存一部分 read 做重叠分析数据量大时内存峰值会很高。解决办法有两个一是降低线程数比如从 8 降到 4内存峰值会明显下降二是用--detect_adapter_for_pe之外的方式比如手动指定接头跳过自动检测阶段。如果你确定文库用的是标准接头手动指定能省下这部分内存。另外输出压缩格式虽然省磁盘但压缩过程也占内存如果内存实在紧张可以先输出未压缩格式后续再压缩。5.2 压缩格式与管道操作的注意事项Fastp 支持直接读写.gz文件但如果你用管道把其他工具的输出接进来要注意格式匹配。比如zcat sample.fastq.gz | fastp -i /dev/stdin ...这种写法Fastp 从标准输入读但输出如果也走标准输出报告文件就没法生成了。更稳妥的做法是用命名管道或者临时文件。如果一定要用管道确保输入输出格式一致并且给报告文件指定明确的路径。还有一个细节Fastp 对.gz的识别是基于文件扩展名的如果你用管道输入它可能不知道数据是压缩的需要加--in1之类的参数明确指定或者用--stdin参数。5.3 UMI 处理中的常见错误UMI 处理最容易出问题的地方是位置和长度不匹配。如果你的 UMI 实际是 10bp但参数写了 8bp切出来的 UMI 就不完整后续去重会出错。反过来如果 UMI 是 8bp 但写了 10bp会把插入片段的前 2bp 也切进去同样有问题。排查方法是先跑一个小样本看输出 read 名称里的 UMI 序列。Fastp 会把 UMI 加到 read 名称里格式是UMI_长度_序列。你可以用zcat sample.clean.fastq.gz | head -4看前几条 read 的名称确认 UMI 长度和内容是否符合预期。如果不对调整--umi_len重新跑。还有一个坑如果 UMI 在 index 里你需要确保 index 读取正确。有些测序平台的 index 是双端 indexread1 和 read2 各有一个 index这时候要用--umi_locper_index并配合--umi_len指定每个 index 的长度。这个配置比较复杂建议先用小样本测试。5.4 过滤过度与过滤不足的平衡这是最考验经验的地方。过滤太狠数据量不够下游分析统计效力下降过滤太松低质量数据混进去比对率低、假阳性高。我的经验是先跑一次默认参数看报告里的各项指标然后根据下游需求调整。如果下游是 SNP calling对错误率敏感可以适当收紧质量过滤-q提到 20 到 25-M提到 25 到 28。如果下游是表达定量对 read 数量更敏感可以放宽到-q 15、-M 20保证数据量。如果下游是组装长度比质量更重要-l可以设低一点比如 25让更多短 read 参与组装。我一般会保留一份默认参数的输出和一份收紧参数的输出对比下游结果看哪个更符合预期。这个对比过程虽然多花时间但能帮你找到最适合自己项目的参数组合。5.5 报告文件路径与权限问题批量脚本里如果报告输出路径没建好Fastp 会直接报错退出。我见过有人脚本里写了-h reports/${sample}.html但忘了mkdir reports结果第一个样本就失败。解决办法是在脚本开头统一建目录并且用-p参数确保目录存在时不报错。权限问题也常见。如果输出目录没有写权限Fastp 会报 permission denied。在共享服务器上跑的时候先确认当前用户对输出目录有写权限。另外如果输入文件是别人创建的读权限也要确认。这些看起来是小事但在批量跑的时候会浪费很多排查时间。6. 把 Fastp 嵌进完整分析流程的几点经验Fastp 很少单独使用它通常是分析流程的第一步。怎么把它和上下游工具衔接好有几个实际经验值得分享。6.1 与比对工具的衔接Fastp 输出的 clean fastq 直接喂给比对工具比如 BWA、Bowtie2、HISAT2。这里要注意 read 名称的一致性。Fastp 默认会保留原始 read 名称但如果启用了 UMI 处理名称会加上 UMI 前缀。比对工具通常不关心名称格式但后续去重工具需要根据 UMI 识别同一分子所以名称格式要统一。另外Fastp 输出的 fastq 质量值编码默认是 Phred33和大多数比对工具兼容。如果你用的工具需要 Phred64需要额外转换但这种情况现在很少见了。6.2 与去重工具的配合如果做了 UMI 处理去重是必须的一步。常用的去重工具有 Picard MarkDuplicates、samtools markdup、UMI-tools 等。这些工具对 UMI 的识别方式不同有的从 read 名称里提取有的需要单独的 UMI 文件。Fastp 把 UMI 放在 read 名称里大多数工具都能识别但格式要匹配。我一般会在 Fastp 之后加一步检查确认 UMI 正确嵌入 read 名称然后再进入比对和去重。如果 UMI 格式不对去重工具可能识别不到导致去重失效。这个检查用zcat | head看几条 read 名称就能完成花不了多少时间但能避免后面的大麻烦。6.3 流程自动化与日志管理在批量项目里我习惯把 Fastp 封装成一个函数或者独立脚本接受样本名和路径作为参数输出固定的目录结构。这样整个流程可以串起来从原始数据到最终结果一条命令跑完。日志管理也很重要。Fastp 的标准输出和标准错误里包含运行信息我一般重定向到每个样本独立的日志文件方便出问题时回溯。日志文件按样本名命名放在统一的 logs 目录下。如果项目很大还可以加时间戳记录每个样本的处理时间方便评估流程效率。6.4 版本控制与参数记录Fastp 的版本更新会带来行为变化比如默认参数调整或者新功能加入。我在项目里会记录使用的 Fastp 版本号用fastp --version获取写进流程文档。这样如果结果有异常可以排查是不是版本差异导致的。参数记录同样重要。我会把每个项目用的完整命令保存成脚本文件和结果一起归档。这样半年后回头看能清楚知道当时用了什么参数为什么这么选。这个习惯在需要复现结果或者排查问题时特别有用。7. 一些让效率翻倍的小技巧最后分享几个我在实际使用中总结的小技巧都是能直接提升效率的。第一个是预检数据量。在跑 Fastp 之前先用zcat sample.fastq.gz | wc -l除以 4 估算 read 数或者用ls -lh看文件大小。这样你能预估处理时间合理安排任务。如果文件特别大可以先抽一部分 read 测试参数确认没问题再跑全量。第二个是复用报告。Fastp 的 JSON 报告里包含所有统计信息你可以写一个脚本把所有样本的 JSON 汇总成一张表用 R 或者 Python 画图。这样不用一个个打开 HTML就能快速看出哪些样本异常。第三个是参数模板化。把常用的参数组合写成几个模板比如严格模式、宽松模式、RNA-seq 模式根据项目类型直接套用。这样既保证一致性又减少每次重新想参数的时间。第四个是定期清理中间文件。Fastp 的输出文件加上报告文件一个样本可能占几个 GB。项目跑完后确认结果没问题及时清理不需要的中间文件避免磁盘爆满。我一般保留 clean fastq 和报告原始数据根据项目要求决定是否保留。第五个是用nohup或screen跑长任务。批量处理可能跑几个小时如果终端断开任务就中断了。用nohup bash run_fastp.sh 或者开一个screen会话能保证任务在后台稳定运行。跑完之后再回来检查日志和结果。这些技巧看起来简单但在实际项目里能省下大量时间和精力。Fastp 本身是个很成熟的工具把参数理解透、把流程组织好它就能成为你分析流程里最可靠的一环。
返回列表