ARTICLE DETAIL

资讯详情

深耕编程入门与网站建设的一线实战洞察。

Fastp实战指南:从参数配置到批量处理,全面掌握fastq质控流程

Fastp实战指南:从参数配置到批量处理,全面掌握fastq质控流程 1. 为什么我最终把质控流程全压在了Fastp上测序数据下机之后第一件事永远是质控。这个环节做得好不好直接决定了后续比对率、变异检出灵敏度、甚至整个项目的成败。我见过太多人拿到fastq文件之后随便跑个FastQC看一眼觉得差不多还行就直接进入比对流程结果后面各种报错、各种异常回头排查半天发现根子还是在原始数据上。Fastp这个工具我从它刚出来没多久就开始用到现在基本上已经成了我处理fastq数据的首选质控方案。原因很简单它把过去需要FastQC Trimmomatic cutadapt 自写脚本才能搞定的事情全部集成到了一个命令里而且速度极快多线程支持好输出报告直观还能自动检测adapter序列。对于日常的WGS、WES、RNA-seq数据Fastp基本能覆盖90%以上的质控需求。这篇文章我打算把自己这些年用Fastp的实战经验完整梳理一遍从最基础的单端/双端处理到UMI处理、polyG修剪、去重、与Shell脚本的配合使用再到那些官方文档里不会写的坑和技巧。适合刚接触fastq质控的新手也适合已经用过Fastp但想进一步榨干它性能的老手。读完你至少能做到拿到一批fastq文件知道该怎么设计质控方案参数怎么调出了问题怎么排查。2. Fastp到底解决了什么问题核心设计思路拆解2.1 传统质控流程的痛点在哪里早些年做质控标准流程是这样的先用FastQC跑一遍原始数据看质量分布、GC含量、adapter污染情况然后根据FastQC的报告手动决定用Trimmomatic还是cutadapt做修剪修剪完再用FastQC跑一遍确认效果如果要做UMI处理还得额外写脚本或者用fgbio、umi_tools这些工具。整个流程下来光是工具之间的格式转换、参数协调就够折腾的。更麻烦的是Trimmomatic的adapter序列需要你自己提供不同建库试剂盒的adapter序列还不一样万一填错了修剪效果大打折扣。cutadapt虽然支持自动检测adapter但速度偏慢而且它只管修剪不管质量过滤和报告生成。Fastp的设计思路就是把这些零散的需求整合起来。它内置了常见的adapter序列库支持自动检测adapter同时完成质量过滤、长度过滤、polyG/polyX修剪、duplication评估、overrepresentation分析最后生成一份HTML报告和JSON格式的统计文件。一个命令搞定全流程这就是它的核心价值。2.2 Fastp的核心功能模块Fastp的功能可以分成几个大块来理解adapter自动检测与修剪。这是Fastp最实用的功能之一。它通过对reads两端进行k-mer分析自动识别adapter序列不需要你手动指定。对于双端数据它还能利用pair overlap信息来辅助判断。实测下来对于标准建库的Illumina数据自动检测的准确率非常高。质量过滤与修剪。Fastp支持按碱基质量进行滑窗修剪类似Trimmomatic的SLIDINGWINDOW也支持按平均质量过滤整条read。默认情况下它会从read两端向中间扫描遇到质量低于阈值的碱基就切掉。这个逻辑比固定长度截断要合理得多。polyG/polyX修剪。这个问题在NextSeq/NovaSeq平台上特别常见因为双色荧光检测的原因信号缺失时会被读成G。Fastp默认开启polyG修剪对这类数据非常友好。polyX修剪则针对其他类型的均聚物污染。UMI处理。Fastp支持在质控阶段直接处理UMI包括UMI的提取、移动到read名称中、以及基于UMI的去重。这个功能对于低起始量建库、ctDNA检测等场景非常关键。重复序列评估。Fastp会统计duplication rate虽然它不做实际的去重操作那是MarkDuplicates或UMI去重的事但这个指标对于判断文库复杂度很有参考价值。过表达序列分析。Fastp会自动检测过表达的序列这对于发现污染、接头二聚体、rRNA残留等问题很有帮助。2.3 为什么选择Fastp而不是其他工具我对比过几个主流方案工具组合速度功能覆盖易用性报告质量FastQC Trimmomatic中等需手动配置adapter一般分开两份报告FastQC cutadapt较慢adapter检测好一般分开两份报告Fastp快全集成一条命令单份HTML报告自写Python脚本慢完全自定义差需自己写Fastp在速度上的优势特别明显。它用C编写多线程效率高处理一个30X的WGS样本约60GB fastq在16线程下大概20-30分钟就能跑完。同样的数据用Trimmomatic时间大概要翻倍。当然Fastp也不是万能的。它不做比对所以不能做基于比对的污染检测它的去重是基于UMI或序列本身的不如Picard MarkDuplicates那么精细。但对于质控这个环节来说Fastp已经足够好了。3. 核心参数详解与实操配置3.1 输入输出参数最基础但最容易出错Fastp的输入输出参数看起来简单但有几个细节不注意就会踩坑。单端数据fastp -i input.fq -o output.fq双端数据fastp -i read1.fq -I read2.fq -o clean1.fq -O clean2.fq注意双端模式下read1和read2的输入输出参数大小写不同-i和-I-o和-O。这个设计一开始让我很不习惯但用多了也就记住了。输出报告fastp -i input.fq -o output.fq -h report.html -j report.json-h指定HTML报告路径-j指定JSON报告路径。JSON报告特别有用方便后续用脚本提取统计指标做批量汇总。注意Fastp默认会输出到标准输出stdout如果你不指定-o结果会直接打印到终端。我第一次用的时候没注意屏幕上刷了一大堆序列吓了一跳。所以务必指定输出文件。还有一个实用参数是--stdout当你需要把Fastp接入管道时很有用fastp -i input.fq --stdout | other_tool3.2 质量过滤参数怎么设才合理Fastp的质量过滤参数主要有这几个-q/--qualified_quality_phred碱基质量阈值默认15。意思是质量值低于15的碱基会被认为是不合格。-u/--unqualified_percent_limit一条read中不合格碱基的比例上限默认40%。超过这个比例整条read被丢弃。-n/--n_base_limit一条read中N碱基的数量上限默认5。-l/--length_required修剪后read的最小长度默认15。--cut_front/--cut_tail/--cut_right从不同方向进行质量修剪。关于-q的设置我的经验是对于Illumina数据Q15作为过滤阈值偏宽松Q20更严格一些。但也不要设太高Q30会导致大量数据被丢弃尤其是read末端。一般WGS用Q15RNA-seq用Q20具体要看数据质量。--cut_right这个参数值得单独说一下。它开启滑窗修剪模式窗口大小由--cut_window_size控制默认4滑窗平均质量低于--cut_mean_quality默认20时从窗口开始位置截断。这个逻辑和Trimmomatic的SLIDINGWINDOW很像但Fastp的实现更快。fastp -i input.fq -o output.fq \ -q 20 -u 30 -n 3 -l 36 \ --cut_right --cut_window_size 4 --cut_mean_quality 20这套参数是我常用的严格模式适合对数据质量要求高的场景。3.3 Adapter修剪参数自动还是手动Fastp默认开启adapter自动检测和修剪。相关参数-a/--adapter_sequence指定read1的adapter序列。--adapter_sequence_r2指定read2的adapter序列。--detect_adapter_for_pe双端模式下开启adapter自动检测。这里有个坑双端模式下adapter自动检测默认是关闭的你需要显式加上--detect_adapter_for_pe。这个设计是因为双端adapter检测需要额外的计算Fastp为了速度默认关掉了。但实际使用中我强烈建议双端数据都加上这个参数。fastp -i r1.fq -I r2.fq -o c1.fq -O c2.fq \ --detect_adapter_for_pe如果你明确知道adapter序列手动指定会更快更准fastp -i r1.fq -I r2.fq -o c1.fq -O c2.fq \ -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCA \ --adapter_sequence_r2 AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT这是标准Illumina TruSeq的adapter序列。不同试剂盒的adapter可能不同用之前最好确认一下。还有一个参数--adapter_fasta可以提供一个FASTA文件包含多个adapter序列Fastp会逐一比对。对于可能混有多种adapter的样本很有用。3.4 polyG和polyX修剪平台特异性处理-g/--trim_poly_g开启polyG修剪默认开启。--poly_g_min_len控制最小长度默认10。-x/--trim_poly_x开启polyX修剪默认关闭。--poly_x_min_len默认10。NextSeq和NovaSeq的数据一定要保持polyG修剪开启。我遇到过好几次用户拿着NovaSeq数据来找我说比对率特别低一看fastqread末端全是G明显是polyG问题。开启修剪后比对率直接上去了。对于其他平台的数据polyX修剪可以按需开启。但要注意polyX修剪可能会误伤真实的polyA尾巴比如RNA-seq中所以RNA-seq数据要谨慎使用。3.5 UMI处理参数低起始量建库的救星UMIUnique Molecular Identifier是一段随机序列在建库时加到每个原始DNA分子上。它的作用是标记每个原始分子这样在后续分析中可以区分真正的生物学重复和PCR扩增产生的重复。Fastp处理UMI的方式很灵活支持多种UMI位置-U/--umi开启UMI处理。--umi_locUMI的位置可选index1、index2、read1、read2、per_index、per_read。--umi_lenUMI长度。--umi_prefixUMI前缀添加到read名称中。--umi_skipUMI后面需要跳过的碱基数。最常见的场景是UMI在read1的前几个碱基fastp -i r1.fq -I r2.fq -o c1.fq -O c2.fq \ -U --umi_locread1 --umi_len8 --umi_prefixUMI这样处理后UMI序列会被移到read名称中格式类似UMI_ACGTACGT:...。后续用umi_tools或fgbio做去重时就能直接从这个名称中提取UMI。如果UMI在index中fastp -i r1.fq -I r2.fq -o c1.fq -O c2.fq \ -U --umi_locindex1 --umi_len8注意Fastp的UMI处理只是把UMI提取出来放到read名称中它本身不做基于UMI的去重。去重需要后续用专门工具完成。这一点很多人会误解。3.6 去重参数什么时候该开Fastp提供了-D/--dedup参数开启基于序列的去重。它的逻辑是如果两条read的序列完全相同或者UMI相同只保留一条。这个功能对于去除PCR重复有一定效果但要注意它不基于比对位置所以不同基因组位置但序列相同的read会被误判为重复。对于高覆盖度数据去重会显著减少数据量。如果后续还要用MarkDuplicates这里去重就重复了。我的建议是如果做了UMI用--dedup基于UMI去重是合理的如果没有UMI一般不在Fastp阶段去重留给后续的MarkDuplicates处理。fastp -i r1.fq -I r2.fq -o c1.fq -O c2.fq \ -U --umi_locread1 --umi_len8 --dedup3.7 其他实用参数-w/--thread线程数默认3。建议设为CPU核心数的1-1.5倍。我一般用16或32。-z/--compression输出压缩级别1-9默认4。级别越高压缩率越好但越慢。4是个不错的平衡点。--overrepresentation_analysis开启过表达序列分析。对于发现污染很有用但会增加运行时间。--correction开启碱基校正。利用双端overlap信息校正错误碱基。对于低质量数据有帮助但也会增加时间。-R/--report_title自定义报告标题。批量处理时很有用方便区分不同样本。4. 完整实操流程从原始数据到干净fastq4.1 场景设定与数据准备假设我们有一批双端测序数据来自NovaSeq平台建库时加了8bp的UMI在read1前端需要做完整的质控处理。数据存放在raw_data/目录下样本名为sample1到sample10。先看一下数据结构ls raw_data/ # sample1_R1.fastq.gz sample1_R2.fastq.gz # sample2_R1.fastq.gz sample2_R2.fastq.gz # ...4.2 单样本Fastp命令设计针对这个场景我设计的Fastp命令如下fastp \ -i raw_data/sample1_R1.fastq.gz \ -I raw_data/sample1_R2.fastq.gz \ -o clean_data/sample1_R1.clean.fastq.gz \ -O clean_data/sample1_R2.clean.fastq.gz \ -w 16 \ -q 20 \ -u 30 \ -n 3 \ -l 36 \ --detect_adapter_for_pe \ --cut_right \ --cut_window_size 4 \ --cut_mean_quality 20 \ -g \ --poly_g_min_len 10 \ -U \ --umi_locread1 \ --umi_len8 \ --umi_prefixUMI \ --dedup \ -h reports/sample1.html \ -j reports/sample1.json \ -R sample1 QC Report \ -z 4逐段解释这个命令的设计逻辑输入输出用-i/-I指定双端输入-o/-O指定双端输出。输出用.gz压缩节省磁盘空间。线程-w 16根据服务器配置调整。Fastp的多线程效率很高16线程基本能跑满。质量过滤-q 20 -u 30 -n 3 -l 36。Q20阈值不合格碱基比例上限30%N碱基上限3最小长度36。这套参数比默认值严格适合对质量要求高的项目。adapter检测--detect_adapter_for_pe双端模式必须显式开启。滑窗修剪--cut_right配合窗口大小4和平均质量20。从5到3扫描遇到低质量窗口就截断。polyG修剪-g --poly_g_min_len 10NovaSeq数据必开。UMI处理-U --umi_locread1 --umi_len8 --umi_prefixUMI从read1前端提取8bp UMI。去重--dedup基于UMI去重。报告-h和-j分别输出HTML和JSON报告-R设置报告标题。压缩-z 4平衡压缩率和速度。4.3 批量处理的Shell脚本单个样本跑通了接下来要批量处理。这里用Shell脚本的for循环来实现#!/bin/bash # 配置 RAW_DIRraw_data CLEAN_DIRclean_data REPORT_DIRreports THREADS16 # 创建输出目录 mkdir -p ${CLEAN_DIR} ${REPORT_DIR} # 获取样本列表 for R1 in ${RAW_DIR}/*_R1.fastq.gz; do # 提取样本名 sample$(basename ${R1} _R1.fastq.gz) R2${RAW_DIR}/${sample}_R2.fastq.gz # 检查R2是否存在 if [ ! -f ${R2} ]; then echo Warning: ${R2} not found, skipping ${sample} continue fi echo Processing ${sample}... fastp \ -i ${R1} \ -I ${R2} \ -o ${CLEAN_DIR}/${sample}_R1.clean.fastq.gz \ -O ${CLEAN_DIR}/${sample}_R2.clean.fastq.gz \ -w ${THREADS} \ -q 20 -u 30 -n 3 -l 36 \ --detect_adapter_for_pe \ --cut_right --cut_window_size 4 --cut_mean_quality 20 \ -g --poly_g_min_len 10 \ -U --umi_locread1 --umi_len8 --umi_prefixUMI \ --dedup \ -h ${REPORT_DIR}/${sample}.html \ -j ${REPORT_DIR}/${sample}.json \ -R ${sample} QC Report \ -z 4 if [ $? -eq 0 ]; then echo ${sample} done. else echo Error: ${sample} failed! fi done echo All samples processed.这个脚本有几个细节值得说明basename ${R1} _R1.fastq.gz用来提取样本名。basename的第二个参数是后缀它会去掉这个后缀。比如raw_data/sample1_R1.fastq.gz会变成sample1。if [ ! -f ${R2} ]检查R2文件是否存在避免因为缺失文件导致脚本中断。if [ $? -eq 0 ]检查上一条命令的退出状态。Fastp成功返回0失败返回非0。这样可以及时发现失败的样本。4.4 并行加速让批量处理快起来上面的脚本是串行执行的10个样本要跑10次。如果服务器核心多可以并行处理。有两种方式方式一用GNU parallells raw_data/*_R1.fastq.gz | \ sed s/_R1.fastq.gz// | \ parallel -j 4 sample{} 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 \ -w 8 -q 20 -u 30 -n 3 -l 36 \ --detect_adapter_for_pe \ --cut_right --cut_window_size 4 --cut_mean_quality 20 \ -g --poly_g_min_len 10 \ -U --umi_locread1 --umi_len8 --umi_prefixUMI \ --dedup \ -h reports/${sample}.html \ -j reports/${sample}.json \ -R ${sample} QC Report -z 4 -j 4表示同时跑4个样本每个样本用8线程。总共32线程。这样比串行快很多。方式二用Shell的后台任务for R1 in raw_data/*_R1.fastq.gz; do sample$(basename ${R1} _R1.fastq.gz) R2raw_data/${sample}_R2.fastq.gz fastp -i ${R1} -I ${R2} \ -o clean_data/${sample}_R1.clean.fastq.gz \ -O clean_data/${sample}_R2.clean.fastq.gz \ -w 8 -q 20 -u 30 -n 3 -l 36 \ --detect_adapter_for_pe \ --cut_right --cut_window_size 4 --cut_mean_quality 20 \ -g --poly_g_min_len 10 \ -U --umi_locread1 --umi_len8 --umi_prefixUMI \ --dedup \ -h reports/${sample}.html \ -j reports/${sample}.json \ -R ${sample} QC Report -z 4 # 控制并发数 while [ $(jobs -r | wc -l) -ge 4 ]; do sleep 1 done done wait echo All done.把命令放到后台执行jobs -r | wc -l统计正在运行的后台任务数超过4个就等待。wait等待所有后台任务完成。注意并行处理时要注意内存和磁盘I/O。Fastp本身内存占用不大但多个实例同时读写磁盘可能会成为瓶颈。如果发现速度上不去可能是磁盘I/O限制了。4.5 结果验证与报告解读跑完之后先看JSON报告里的关键指标# 提取关键统计 cat reports/sample1.json | python3 -c import json, sys data json.load(sys.stdin) summary data[summary] before summary[before_filtering] after summary[after_filtering] print(f\Total reads before: {before[total_reads]}\) print(f\Total reads after: {after[total_reads]}\) print(f\Pass rate: {after[total_reads]/before[total_reads]*100:.2f}%\) print(f\Q30 before: {before[q30_rate]*100:.2f}%\) print(f\Q30 after: {after[q30_rate]*100:.2f}%\) print(f\GC content: {after[gc_content]*100:.2f}%\) 一般来说合格的质控结果应该满足Pass rate在80%以上严格模式可能70%左右Q30 after在90%以上GC content与物种预期相符Duplication rate在合理范围WGS一般20%RNA-seq可能更高如果Pass rate过低要检查是不是质量阈值设得太严或者原始数据本身质量差。如果Q30 after没有明显提升说明修剪效果不好可能需要调整参数。HTML报告更直观用浏览器打开就能看到各种图表质量分布、碱基组成、adapter含量、duplication等。我一般会重点看几个图Quality curve修剪前后的质量分布对比Base composition碱基组成是否正常有没有异常偏移Adapter contentadapter是否被有效去除Duplication重复率是否正常5. 常见问题与排查技巧实录5.1 Fastp运行报错怎么办报错一Failed to open file最常见的原因是路径写错了或者文件不存在。检查ls -lh raw_data/sample1_R1.fastq.gz如果文件存在但还是报错可能是权限问题chmod r raw_data/sample1_R1.fastq.gz报错二std::bad_alloc内存不足。Fastp处理大文件时需要一定内存尤其是开启overrepresentation分析时。解决方法减少线程数每个线程会占用一定内存关闭overrepresentation分析增加服务器内存报错三Segmentation fault段错误通常是输入文件损坏。检查fastq文件完整性gzip -t raw_data/sample1_R1.fastq.gz如果gzip测试报错说明文件损坏需要重新下载或从备份恢复。5.2 质控效果不理想的排查思路问题adapter去除不干净先确认adapter序列是否正确。如果自动检测效果不好手动指定fastp -i r1.fq -I r2.fq -o c1.fq -O c2.fq \ -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCA \ --adapter_sequence_r2 AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT如果还是不行可能是adapter发生了突变或者有多个adapter。用--adapter_fasta提供多个候选序列。问题polyG修剪过度如果发现read被切得太短可能是polyG阈值设得太低。默认--poly_g_min_len 10可以提高到15或20fastp -i r1.fq -I r2.fq -o c1.fq -O c2.fq \ -g --poly_g_min_len 15问题UMI提取不正确检查UMI位置和长度是否与建库方案一致。如果不确定可以先不开启UMI处理跑一遍Fastp看看read前几个碱基的组成。如果前8bp的碱基组成是随机的每种碱基约25%那很可能就是UMI。# 查看前10bp的碱基组成 zcat raw_data/sample1_R1.fastq.gz | head -1000 | \ awk NR%42 {print substr($0,1,10)} | \ sort | uniq -c | sort -rn | head -205.3 性能优化技巧技巧一用管道避免中间文件如果Fastp后面还有其他处理步骤可以用管道直接传递fastp -i input.fq --stdout -w 8 | \ bwa mem -t 8 reference.fa - | \ samtools sort - 8 -o output.bam这样避免了写中间fastq文件节省磁盘I/O。技巧二合理设置压缩级别-z 4是默认值平衡了速度和压缩率。如果磁盘空间紧张可以设-z 6或-z 9但会慢一些。如果追求速度设-z 1或-z 2。技巧三用--stdin从管道读取Fastp支持从标准输入读取zcat input.fq.gz | fastp --stdin -o output.fq这在处理流式数据时很有用。5.4 常见问题速查表问题现象可能原因解决方法Pass rate过低质量阈值太严降低-q或-uQ30 after无提升修剪参数不当调整--cut_right参数adapter残留检测未开启或序列不对加--detect_adapter_for_pe或手动指定read过短polyG/polyX修剪过度提高--poly_g_min_lenUMI未提取位置或长度不对检查--umi_loc和--umi_len运行速度慢线程数不足或I/O瓶颈增加-w检查磁盘内存不足线程过多或开启overrepresentation减少线程关闭overrepresentation输出文件为空输入文件损坏或参数错误检查输入文件查看日志5.5 几个我踩过的坑坑一双端adapter检测默认关闭这个前面提过但值得再强调。我第一次用Fastp处理双端数据时没加--detect_adapter_for_pe结果adapter残留严重。后来看文档才发现这个默认行为。现在我的脚本里这个参数是必加的。坑二UMI处理不等于UMI去重Fastp的-U只是把UMI提取到read名称中--dedup才是去重。而且--dedup在没有UMI时是基于序列去重有UMI时是基于UMI去重。这两个参数要配合使用。坑三polyG修剪对非NextSeq数据的影响有一次处理MiSeq数据习惯性开了polyG修剪结果发现一些真实的polyG区域被误切了。后来查资料才知道polyG问题主要是双色荧光平台的特性MiSeq是四色荧光不存在这个问题。所以polyG修剪要按平台决定是否开启。坑四JSON报告中的百分比是小数Fastp的JSON报告中q30_rate、gc_content这些字段是0-1之间的小数不是百分比。我第一次用脚本提取时忘了乘100结果报告里Q30是0.95%闹了笑话。坑五输出文件名不要和输入相同Fastp会先读输入再写输出但如果输入输出文件名相同可能会出问题。虽然Fastp有保护机制但最好还是用不同的文件名。6. 进阶技巧让Fastp发挥更大价值6.1 与umi_tools的配合使用Fastp处理完UMI后read名称中会带有UMI信息。接下来可以用umi_tools做基于UMI的去重和定量# Fastp处理 fastp -i r1.fq -I r2.fq -o c1.fq -O c2.fq \ -U --umi_locread1 --umi_len8 --umi_prefixUMI # 比对 bwa mem -t 16 reference.fa c1.fq c2.fq | \ samtools sort - 16 -o aligned.bam # 提取UMI并去重 umi_tools extract --stdinc1.fq --stdoutc1.umi.fq \ --read2-inc2.fq --read2-outc2.umi.fq \ --bc-patternNNNNNNNN # 或者从BAM中提取 umi_tools dedup -I aligned.bam -S dedup.bam \ --extract-umi-methodread_name \ --umi-separator:注意Fastp的UMI前缀格式是UMI_umi_tools默认的分隔符是_可能需要调整--umi-separator参数。6.2 用Shell脚本做质控报告汇总批量处理完后通常需要把所有样本的统计指标汇总成一张表。用Shell Python可以轻松搞定#!/bin/bash echo -e Sample\tTotal_Reads\tPass_Rate\tQ30_Before\tQ30_After\tGC\tDup_Rate qc_summary.tsv for json in reports/*.json; do sample$(basename ${json} .json) python3 -c import json with open(${json}) as f: data json.load(f) s data[summary] b s[before_filtering] a s[after_filtering] d data.get(duplication, {}) print(f\${sample}\t{b[total_reads]}\t{a[total_reads]/b[total_reads]*100:.2f}%\t{b[q30_rate]*100:.2f}%\t{a[q30_rate]*100:.2f}%\t{a[gc_content]*100:.2f}%\t{d.get(rate, 0)*100:.2f}%\) qc_summary.tsv done echo Summary written to qc_summary.tsv这个脚本会生成一个TSV文件可以直接用Excel打开或者导入R做可视化。6.3 处理特殊类型数据的参数调整RNA-seq数据不要开polyX修剪会误伤polyA尾巴--cut_right可以保留但窗口质量阈值可以放宽到15考虑开启--correction利用双端overlap校正错误fastp -i r1.fq -I r2.fq -o c1.fq -O c2.fq \ --detect_adapter_for_pe \ --cut_right --cut_window_size 4 --cut_mean_quality 15 \ -g --poly_g_min_len 10 \ --correction \ -q 15 -u 40 -l 36扩增子数据长度过滤要严格因为扩增子长度是已知的adapter修剪要彻底因为扩增子数据adapter污染通常较严重考虑开启--dedup扩增子数据重复率通常很高fastp -i r1.fq -I r2.fq -o c1.fq -O c2.fq \ --detect_adapter_for_pe \ -l 200 -q 20 -u 20 \ --dedup \ -g低起始量数据UMI处理必开去重必开质量阈值可以适当放宽因为低起始量数据本身质量可能就一般fastp -i r1.fq -I r2.fq -o c1.fq -O c2.fq \ --detect_adapter_for_pe \ -U --umi_locread1 --umi_len8 --umi_prefixUMI \ --dedup \ -q 15 -u 40 -l 30 \ -g6.4 与Nextflow/Snakemake的集成如果项目规模大建议把Fastp集成到工作流管理工具中。以Nextflow为例process fastp { tag $sample publishDir results/clean, mode: copy input: tuple val(sample), path(r1), path(r2) output: tuple val(sample), path(${sample}_R1.clean.fastq.gz), path(${sample}_R2.clean.fastq.gz) path(${sample}.html) path(${sample}.json) script: fastp \\ -i ${r1} -I ${r2} \\ -o ${sample}_R1.clean.fastq.gz \\ -O ${sample}_R2.clean.fastq.gz \\ -w ${task.cpus} \\ -q 20 -u 30 -n 3 -l 36 \\ --detect_adapter_for_pe \\ --cut_right --cut_window_size 4 --cut_mean_quality 20 \\ -g --poly_g_min_len 10 \\ -U --umi_locread1 --umi_len8 --umi_prefixUMI \\ --dedup \\ -h ${sample}.html -j ${sample}.json \\ -R ${sample} QC Report -z 4 }这样可以利用Nextflow的并行调度、错误重试、断点续跑等功能比裸写Shell脚本更可靠。6.5 质控后的数据管理建议质控完成后建议保留以下文件干净的fastq文件压缩Fastp的HTML报告Fastp的JSON报告运行日志如果有使用的命令或脚本目录结构建议project/ ├── raw_data/ # 原始数据只读 ├── clean_data/ # 质控后数据 ├── reports/ # 质控报告 │ ├── html/ │ └── json/ ├── scripts/ # 处理脚本 ├── logs/ # 运行日志 └── qc_summary.tsv # 汇总统计原始数据一定要保留不要为了省空间删掉。万一后续发现质控参数需要调整还能重新跑。我一般会在项目结束后把原始数据归档到冷存储但至少保留到项目发表。7. 一些个人体会Fastp这个工具我用了快五年从最早的0.20版本到现在功能越来越完善。它最大的价值在于把质控流程标准化了——以前每个人写的质控脚本都不一样现在一个Fastp命令就能覆盖大部分需求结果可复现报告可比较。但工具再好也只是工具。质控参数怎么设还是要根据具体的数据类型、建库方案、下游分析需求来决定。我见过有人直接抄别人的参数结果数据被切得太狠后面分析灵敏度下降。也见过有人参数太宽松adapter没去干净导致假阳性。我的建议是新项目开始时先拿一两个样本做参数测试看看不同参数下的Pass rate、Q30、adapter去除效果找到最适合这个项目的参数组合再批量处理。这个前期投入是值得的能避免后面返工。另外Fastp的报告一定要认真看。很多人跑完Fastp就直接进入下一步报告看都不看。其实报告里有很多有价值的信息过表达序列可能提示污染duplication rate可能提示文库复杂度问题adapter含量可能提示建库质量。这些信息对于判断数据是否可用非常关键。最后说一个实际工作中的小习惯我会把每次Fastp运行的完整命令记录在一个commands.log文件里包括日期、样本名、参数。这样半年后回头看还能知道当时是怎么处理的。这个习惯帮我省了很多次当时到底怎么跑的的困惑。
返回列表