新一批测序数据下机了“先跑一下fastqc”这句话几乎是我每次拿到原始fastq时脑子里蹦出来的第一反应。FastQC是生信圈子里最常用的测序数据质控软件没有之一。我见过太多入门的朋友下载完fastq.gz就直接去比对、定量折腾一整天后结果明显不对劲回头一查才发现原始数据里全是接头序列。这篇文章就把第七课的内容完整讲一遍——从FASTQ格式怎么读到FastQC怎么装、怎么跑、报告怎么看最后再讲质控不过关的数据怎么处理争取让零基础的朋友也能把这一步走扎实。1. 测序数据到手后为什么必须先把“质控”排在所有分析之前1.1 下机fastq里常见的“坑”低质量碱基、接头残留、重复序列很多第一次接触测序数据的朋友会有一个错觉测序仪出来的数据应该是“干净”的毕竟一台仪器几十万上百万出来的数据总不至于太差吧。但实际情况远不是这样。测序仪下机的原始fastq里常见的“坑”至少有这几类第一是低质量碱基。测序过程中荧光信号衰减、焦磷酸测序中的化学噪音、芯片上的局部气泡都会让某一段read的碱基质量掉得很厉害尤其是读长末端质量曲线几乎必然往下走。第二是接头残留。建库的时候DNA片段两端要接上已知序列的接头才能被测序仪识别。如果片段长度比测序读长还短测序仪就会一路读过去把接头序列也读进来。这些接头序列不是目标基因组的一部分属于典型的“污染物”。第三是重复序列抬高。PCR扩增过程中某些片段被过度扩增导致文库中大量read完全一样。这种重复会严重干扰后续的表达定量和变异检测。第四是N碱基过多。测序仪在某些位置实在判断不出是A、T、C还是G就会输出一个N。少量N可以接受N多了说明这部分测序彻底失败。这些“坑”如果不提前发现后面每一步分析都会被污染数据带偏。1.2 不做质控直接比对下游结果会受到什么影响有人会说不就是多几个低质量碱基嘛比对软件比如BWA、STAR本身也会有错配容忍度何必较真如果你只是跑通一个流程看看结果长什么样那确实不较真也能跑。但只要你做的是正经科研项目或者要拿数据发文章质控这步跳过后果会在下游被不断放大接头序列比对不上参考基因组比对率会莫名下降让你误以为物种匹配出了问题。低质量碱基会引起大量错配变异检测时出现一堆假阳性SNP。重复序列会严重高估某些区域的覆盖度影响拷贝数变异和peak calling。接头残留如果出现在RNA-seq数据里还会干扰转录本定量的准确性尤其是短片段。说白了质控不是“洁癖”而是在源头堵住系统性误差。数据一旦送进去跑完整个流程再回头排查就非常被动。1.3 FastQC回答的三个核心问题能不能用、哪里有问题、要不要修FastQC这类软件存在的意义就是快速、自动化地回答三个问题。能不能用看整体质量分布和基本统计如果数据整体质量很差就要重新考虑实验或测序策略。哪里有问题通过十多个模块逐项扫描定位具体问题类型比如是接头污染、重复率高还是某个测序lane出现系统性偏差。要不要修FastQC本身不做修剪但它的结果会告诉你哪些问题需要交给Trimmomatic、cutadapt这类工具去处理哪些问题修不了只能重测。所以FastQC更像是“体检报告”不是“治病药方”。你要学会读报告然后决定下一步怎么办。2. 拆开FASTQ文件读懂四行记录与Phred质量值体系2.1 一条read的四行结构说FastQC之前得先把FASTQ格式讲明白因为FastQC的所有判断都是基于FASTQ的四个信息层。FASTQ格式每个read固定占四行举个实际例子HISEQ:289:C5L7TACXX:6:1101:1234:5678 1:N:0:ATCACG GATTACAGATTACAGATTACAGATTACAGATTACAGATTACAG IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII第一行以开头是read的名称和注释信息第二行是碱基序列本身第三行以开头通常只写一个有时会重复read名称第四行是质量值字符串长度必须与第二行碱基序列完全一致每个字符对应一个碱基的质量编码。很多刚入门的朋友会忽略第二行和第四行长度必须一致这个细节。如果长度对不上说明文件传输损坏或者软件处理出错FastQC一般会直接报错或者警告。2.2 质量值与错误概率的换算逻辑第四行那一串看似乱码的字符其实藏着一个很优雅的数学换算。质量值Q与碱基错误概率P之间的关系是Q -10 × log10(P)所以Q20对应的错误概率是1%Q30对应0.1%Q40对应0.01%。数值越高这个碱基越可靠。我们常说的“测序质量达到Q30”意思是单碱基错误率不超过千分之一。RNA-seq、WGS这类项目通常会要求大部分碱基达到Q30以上。第四行中每个ASCII字符都代表一个数值换算规则是ASCII码值减去一个偏移量就得到Phred质量分数。比如字符‘I’的ASCII码是73如果偏移量是33那质量值就是40。2.3 Phred33还是Phred64编码体系判断方法这里藏着一个老生常谈的坑——不同版本的Illumina测序仪质量值编码体系不一样。Phred33编码ASCII码值 Q 33范围约Q0到Q41对应字符33到74。现在Illumina 1.8版本下机数据基本都是这种。Phred64编码ASCII码值 Q 64范围约Q0到Q62。老版本Illumina 1.3到1.5时代的格式。如果分析软件误判了编码体系质量值会被整体算错导致质控结果完全失真。怎么快速判断最简单的方法就是看质量字符串里出现的字符范围。如果看到‘I’、“H”、“G”这种大写字母密集出现同时还有空格或“#”、“$”这些符号大概率是Phred33如果看到大量小写字母或者、A、B等字符可能要怀疑是Phred64。其实直接看FastQC的Basic Statistics模块它会自动判断并显示Encoding类型这是最省事的办法。2.4 双端文件命名与读长确认拿到fastq数据时还要先确认文件是不是双端测序命名是否规范。常见的双端文件命名是sample_1.fastq.gz sample_2.fastq.gz或者sample_R1.fastq.gz sample_R2.fastq.gzR1和R2的read数量必须一致顺序也要一一对应因为后续比对软件会把它们当作同一DNA片段的两端来处理。如果两个文件的行数不一致或者名字对不上比对的paired-end信息就会错乱。另外下机数据通常还带一个SampleSheet或者测序说明文件里面写着读长是PE150还是PE100。这个信息在后续做读长分布判断时很关键因为实际读长如果明显短于预期就要考虑是不是测序提前终止。3. FastQC的安装、运行与批量操作3.1 安装FastQCconda、官方包与在线版FastQC的安装方式主要有三种我逐个说一下。最推荐的是用conda安装一行命令解决依赖问题conda install -c bioconda fastqc这个方式的好处是conda会自动处理Java依赖。FastQC本质上是一个Java程序老版本要求Java 8新版本甚至要求Java 11以上。如果你直接用官方包而系统里Java版本不对运行时会直接报错conda装就省心很多。第二种是去官网下载Unix版本解压后用目录里的fastqc脚本运行wget https://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.12.1.zip unzip fastqc_v0.12.1.zip cd FastQC chmod x fastqc ./fastqc --version第三种是使用Galaxy在线版适合完全不想装软件的新手直接在网页上传fastq文件点几下就能出报告。但缺点是数据要上传到服务器如果数据量大有隐私要求本地安装更稳妥。我个人的建议是本地能用conda就用conda学习阶段折腾一次环境配置后面所有工具都能跟着受益。3.2 单样本质控命令与常用参数安装完成后最基础的用法就是fastqc sample_R1.fastq.gz默认情况下FastQC会在当前目录生成两个文件sample_R1_fastqc.html和sample_R1_fastqc.zip。实际项目中我更常用带参数的命令fastqc sample_R1.fastq.gz sample_R2.fastq.gz -o fastqc_results -t 4这里解释一下参数含义-o指定输出目录先创建好目录再运行否则FastQC可能直接报错。-t指定线程数一次分析多个文件时可以并行加速。--noextract可以防止自动解压zip包如果你想保留zip里的原始数据文件建议加上。-f fastq显式指定输入格式虽然大部分时候FastQC能自动识别但写清楚更保险。一条命令可以同时传入多个fastq文件。如果你的文件名有规律甚至可以偷懒写fastqc *.fastq.gz -o fastqc_results -t 8通配符会让shell自动展开所有匹配文件FastQC逐个分析一次性出所有报告。3.3 批量处理对一整个目录的fastq.gz跑fastqc真实项目里很少只有一个样本经常是几十上百个样本。这时候不要手动一条条敲命令写个for循环把整个目录扫一遍更靠谱。for f in *.fastq.gz; do fastqc $f -o fastqc_results -t 4 done注意变量名加引号防止文件名带空格时出错。如果样本数量特别多-t线程数不要超过CPU核心数否则反而会拖慢速度。如果追求更规范的批处理还可以这样写把R1和R2放进同一轮循环for f in *_R1.fastq.gz; do r2${f/_R1/_R2} fastqc $f $r2 -o fastqc_results -t 4 done这段脚本会同时处理同一对双端文件最终每个样本生成R1和R2两个html报告。批处理跑完之后建议把所有html文件放到同一个目录里统一打开配合后面的MultiQC汇总会更直观。3.4 输出文件有哪些分别怎么用每个样本分析完会生成两个文件.html可视化报告浏览器直接打开就能逐模块查看。.zip压缩包里面包含fastqc_data.txt、fastqc_report.html、summary.txt等文件。如果你要做自动化统计比如提取每个样本的平均质量值直接解析fastqc_data.txt比读html方便得多。summary.txt里每一行对应一个模块的状态PASS表示通过、WARN表示警告、FAIL表示失败。这个文件非常适合批量脚本汇总我经常用一条awk把它拼成表格for f in fastqc_results/*_fastqc.zip; do unzip -p $f */summary.txt | awk -v name$f {print name, $1, $2} done这一步能把所有样本的各模块状态快速拉出来一眼看出哪几个样本需要重点处理。4. 报告硬指标解读质量值、碱基分布、读长与N含量4.1 Basic Statistics与Per base sequence quality打开html报告最上面的是Basic Statistics模块主要包括Filename文件名File type是常规碱基还是颜色空间编码Encoding质量编码体系比如Sanger / Illumina 1.9看到这个说明是Phred33Total Sequences总read数Sequence length读长比如150 bp%GC整体GC含量这个模块基本没有绿色红色之分主要用来确认样本信息对不对。如果你预期PE150结果多了一个PE300或长度不一致说明样品命名或数据整理出了问题。紧接着是Per base sequence quality这是最常看的模块。它展示每个碱基位置的质量值分布图形中间的红线是中位数黄色盒子是四分位距上下须是10%和90%分位蓝线是平均质量。判断经验很简单如果整个图大部分位置都在绿色区域Q28以上就是一条好曲线如果末端明显掉到Q20以下说明读长末端质量衰减严重如果一开始就很低那整个测序质量都堪忧。FastQC的判定规则大致是如果任何位置的中位数质量低于Q30会警告低于Q20会判失败。注意这里的“失败”只是提示你数据存在质量问题不代表全盘放弃后面修剪仍可能救回来。4.2 Per sequence quality scores与Per tile sequence qualityPer sequence quality scores这个模块展示的是每条read平均质量的分布。横轴是平均质量值纵轴是read数量。正常情况下这个图应该是一个集中在高分数段Q30到Q40的单峰。如果出现一个在低质量区域的拖尾峰说明有一大批read整体质量很差可能是文库制备或测序lane出了问题。FastQC规则中如果最差质量峰低于Q20会失败低于Q30会警告。Per tile sequence quality是一个热图横轴是碱基位置纵轴是测序芯片上的不同tile小区域。这个模块主要看测序芯片不同区域是否存在系统性质量偏差比如某个tile整列颜色异常偏蓝或偏红就说明芯片该区域有气泡或者信号异常。常规小项目里这个模块经常显示一条淡蓝色问题不大如果出现一整块很深的色块需要检查测序服务商是否在这条lane上有局部故障。4.3 Per base sequence content与Per sequence GC contentPer base sequence content展示每个碱基位置上A、T、C、G四种碱基的占比。理论上随机打断的基因组文库中每个位置四种碱基比例应该差不多尤其在read的起始位置不应有明显偏向。如果看到前10-15bp出现强烈的AT或GC分离说明有接头或引物残留导致的偏向性扩增这是常见的建库偏差信号。FastQC的规则是任一位置碱基比例偏差超过10%会警告超过20%会失败。但这里我要提醒一句很多RNA-seq数据或扩增子数据碱基组成本身就不均匀这类模块的警告不一定是数据质量问题要结合实验类型判断。Per sequence GC content展示测序reads的GC含量分布同时叠加一个理论正态分布曲线。如果实测曲线明显偏离理论曲线、出现双峰或偏峰可能意味着文库混入了其他物种DNA、存在PCR偏好性或者样本本身GC含量异常。比如双峰中一个峰在40%另一个在60%大概率是污染或者两个不同来源的DNA混在一起。4.4 碱基N含量与Sequence Length DistributionPer base N content显示每个碱基位置N的比例。N代表测序仪无法判断的碱基通常是因为信号弱或位置重叠。如果N含量在某个位置突然升高比如接近100%往往提示读长末端测序彻底失效。FastQC的规则是N比例超过5%警告超过20%失败。实际项目中正常高质量数据这个模块基本是一条贴近0的直线。如果出现明显的N峰或者N比例整体偏高要考虑是否终止分析或者对末端进行修剪。Sequence Length Distribution显示reads长度的分布。如果是固定读长的Illumina测序通常会看到一个尖锐的单峰比如150 bp。如果长度分散可能是fastq文件混入了不同批次的reads或者测序过程中发生接头连接不均一。FastQC对长度分布的判定比较严格只要存在长度不一致就会警告如果长度低于35 bp会更警惕。对双端数据来说还要注意R1和R2读长是否一致不一致会影响后续比对软件的参数设置。5. 报告脏数据解读重复、过表达序列、接头与Kmer5.1 Sequence Duplication Levels重复多少算异常这个模块展示的是序列重复程度的分布。横轴是重复次数纵轴是reads数量。需要特别说明的是对于WGS全基因组测序一定程度的重复是正常的因为基因组本身包含重复区域高深度测序时即使没有PCR重复随机采样也会产生部分完全相同的reads。但如果重复水平超过一定比例就要警惕了。FastQC的判定规则大致是如果非唯一序列占总序列比例超过20%会警告超过50%会失败。实际经验中WGS数据如果重复率超过30%-40%通常说明文库复杂度偏低可能是起始DNA量不足或者PCR循环数过多。还有一个容易混淆的问题RNA-seq数据通常重复率天然偏高因为高表达基因会产生大量完全相同或几乎相同的reads。这时候看到Duplication模块亮黄灯先别急着下结论可以结合Overrepresented sequences模块一起看。5.2 Overrepresented sequences先查污染再查rRNA这个模块会列出在数据中占比非常高的序列。FastQC认为如果某条序列占总reads数的0.1%以上就算过表达序列会单独列出来并尝试标注它属于什么类型比如是否匹配接头、是否匹配rRNA、是否匹配已知物种序列。看到这个模块报红点第一反应不是害怕而是分类处理如果是接头序列说明接头去除不彻底需要用cutadapt或Trimmomatic再切。如果是rRNA序列说明建库时rRNA去除效率不高这在RNA-seq里也算常见问题但不影响mRNA分析太多。如果是大量未知来源的短序列就要怀疑样本间污染必要时需要做物种来源鉴定。有一点值得注意FastQC只能识别它内置数据库里的常见污染序列未知的污染序列它不一定能标注出来。所以看到“No Hit”的过表达序列时不能直接忽略可以手动把这条序列拿去比对一下。5.3 Adapter Content接头残留如何判断Adapter Content模块展示从哪个碱基位置开始出现接头序列以及累积比例。一张健康的图应该是整条曲线都接近0。如果曲线从某个位置开始突然上升说明reads在这个位置之后开始读到接头。曲线越早上升、最终比例越高问题越严重。FastQC的判定规则大致是接头含量累计超过5%会警告超过10%会失败。实际中最常见的接头残留情况是读长大于插入片段尤其是插入片段比较短的文库比如small RNA-seq或者FFPE样本接头污染会比较明显。处理办法也很成熟Trimmomatic或cutadapt去接头一般用ILLUMINACLIP参数同时做两件事——先识别并切除接头再顺便修掉与接头相邻的低质量碱基。5.4 Kmer Content与另一种观察维度Kmer Content模块在某些FastQC版本中默认不展示但这个模块藏着不少信息。它统计了在特定位置出现频率异常高的短序列通常是4-8bp的kmer并给出富集位置。如果一个kmer只在read开头富集可能是引物或接头部分残留如果它在read中间持续富集可能是某种重复元件如果kmer富集但找不到对应已知序列要怀疑是否出现了未知污染。不过对刚入门的朋友Kmer模块可以先放一放不用一上来就深挖。先把前几个模块读明白日常项目真实需求已经覆盖了90%以上。5.5 12个模块判定规则汇总表为了方便查看和归档我把FastQC报告里常用的模块、状态含义和最常见的应对思路整理成一张表模块警告常见原因失败常见原因常见处理方式Basic Statistics无编码体系无法识别确认测序格式Per base sequence quality末端质量偏低整体质量崩溃末端修剪严重则重测Per tile sequence quality个别tile偏差多个tile系统性偏差联系测序服务商Per sequence quality scores质量峰整体下移大量reads质量极差低质量reads过滤Per base sequence content前几个碱基偏向整段碱基组成失衡检查引物/接头建库偏差Per sequence GC content偏离理论分布出现多峰污染排查检查样本Per base N contentN比例5%N比例20%修整末端或重测Sequence Length Distribution长度不一致长度过短确认读长设置Sequence Duplication Levels重复率偏高重复率极高降低PCR循环数Overrepresented sequences单序列占比0.1%多种污染序列去接头、去rRNA、查污染Adapter Content接头5%接头10%接头修剪Kmer Content特定kmer异常多种kmer异常结合Overrepresented判断这张表是我平时自己看报告时的速查表建议你也存一份尤其是做多个样本对比时很有用。6. 质控不合格怎么办修剪参数、复检流程与我的判据经验6.1 先诊断再处理不同失败类型的应对策略FastQC报告亮红不代表这组数据就废了。首先要诊断是哪种类型的失败。如果是质量值低比如Per base sequence quality末端掉得很厉害可以修剪末端低质量碱基比如SLIDINGWINDOW:4:15这种参数会把滑动窗口内平均质量低于Q15的片段裁掉。如果是接头残留直接切除接头序列可以选择Trimmomatic的ILLUMINACLIP模块或者更轻量的cutadapt。如果是重复率过高修剪解决不了根本问题因为重复是建库阶段引入的只能从实验端想办法但还是要看重复reads是否影响你的分析目标。比如RNA-seq表达定量时有独特的处理方式而WGS做SNP检测时对重复更敏感。如果是N含量过高说明测序信号丢失严重轻度N可以修剪重度N只能重测。所以不要一看到失败就删数据先看失败模块和失败位置再决定“能救”还是“不能救”。6.2 Trimmomatic双端修剪标准流程Trimmomatic是目前最常用的质控修剪工具之一双端数据处理命令我给你一个可以直接改的模板trimmomatic PE \ sample_R1.fastq.gz sample_R2.fastq.gz \ sample_R1_trim.fastq.gz sample_R1_unpaired.fastq.gz \ sample_R2_trim.fastq.gz sample_R2_unpaired.fastq.gz \ ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \ LEADING:3 \ TRAILING:3 \ SLIDINGWINDOW:4:15 \ MINLEN:36逐个参数解释PE表示双端模式输入R1和R2输出四个文件其中*_unpaired是修剪后另一端被淘汰的孤儿reads一般不用。ILLUMINACLIP:TruSeq3-PE.fa:2:30:10是去接头参数冒号分隔的3个数字依次是允许的最大错配数、配对分数阈值、单端分数阈值。实际调优时如果接头残留严重可以适当把阈值调低。LEADING:3和TRAILING:3分别切除5端和3端质量值低于3的碱基。SLIDINGWINDOW:4:15是滑动窗口修剪每4个碱基为一组平均质量低于15就切掉。MINLEN:36是修剪后read长度低于36bp就丢弃。注意TruSeq3-PE.fa是Trimmomatic自带的接头序列文件如果你的建库试剂盒不是Illumina TruSeq要去厂家说明书里找对应接头序列或者手动写成fasta格式加进去。6.3 修剪后复检前后对比怎么看修剪完成后一定要再做一次FastQC看修剪是否到位。复检时重点看三个模块Adapter Content曲线是否降到接近0。如果还是很高说明接头参数没配对或者接头序列文件不对。Per base sequence quality末端是否明显回升。修剪后末端虽然变短但保留下来的碱基质量应该更稳定。Sequence Length Distribution修剪会导致reads长度出现一定波动这是正常的。但如果长度范围过宽要注意MINLEN需要调高。此外我通常还会把质控前后两个报告放到MultiQC里对比multiqc fastqc_before fastqc_after -o multiqc_reportMultiQC会自动扫描目录下所有fastqc报告生成一个汇总交互页面。这样多个样本的状态就能在一张页面上看全不用一个个打开html。6.4 几个判断数据是否可用的个人经验最后分享几条我这几年的实际判断经验供你参考。第一条单看一个模块的FAIL不能下结论要交叉验证。比如Adapter Content FAIL同时Per base sequence content前几个碱基出现偏差那基本可以确定是接头问题但如果Adapter Content FAIL而其他模块全绿可能只是少数reads带接头影响可控。第二条RNA-seq数据要允许适度“不完美”。因为在转录组数据里GC含量不均匀、高表达基因重复高都是生物学特征不是纯技术问题。不要机械地追求12个模块全绿。第三条平均质量Q30以上末端质量掉到Q20附近这在长读长测序中非常常见。只要修剪掉末端后有效数据比例还在70%-80%以上后续分析完全没问题。第四条如果质控后有效数据量骤降比如PE150修剪完只剩下PE100以下那要重新评估测序方案。与其靠修剪硬撑不如重新加大测序深度。培养一个习惯每批数据正式分析前跑一次fastqc修剪完再跑一次中间顺手把MultiQC报告归档。别小看这一步到了写文章或补实验的时候这些质控截图和报告就是最能说明数据质量的证据。生信这个领域很多坑不是靠聪明避开的是靠流程规范一个个提前堵死的。