第3節(jié) 數(shù)據(jù)質(zhì)控
一无切、正式流程的搭建愧杯，整個(gè)完整的流程分為以下6部分：
① 原始測(cè)序數(shù)據(jù) fastq 的質(zhì)控 QC；
② read比對(duì)的圆，排序和去除重復(fù)序列；
③ Indel區(qū)域重（“重新”的“重”）比對(duì)半火；
④ 堿基質(zhì)量值重校正BQSR越妈；
⑤ 變異檢測(cè)；
⑥ 變異結(jié)果質(zhì)控和過(guò)濾VQSR钮糖；
注意：
① 當(dāng)前以 Illumina 為代表的新一代測(cè)序 (NGS) 平臺(tái)梅掠，主要采用邊合成邊測(cè)序（SBS, Sequencing by Synthesis）的技術(shù)。該方法通過(guò)化學(xué)反應(yīng)完成堿基的合成店归，使 DNA 鏈從 5' 端逐步延伸至 3' 端阎抒。然而，隨著合成鏈長(zhǎng)度的增加消痛，DNA 聚合酶的效率會(huì)逐漸降低且叁，其特異性也會(huì)變差，這就導(dǎo)致了堿基合成的錯(cuò)誤率在鏈延伸后期逐步升高秩伞。這種現(xiàn)象是目前 NGS 測(cè)序讀長(zhǎng)普遍較短的主要原因之一逞带。此外，在測(cè)序開(kāi)始階段纱新，由于化學(xué)反應(yīng)尚未穩(wěn)定展氓，可能出現(xiàn)質(zhì)量值波動(dòng)的情況。但這種波動(dòng)通常局限于高質(zhì)量值區(qū)域脸爱，影響相對(duì)較小遇汞。
② 而測(cè)序數(shù)據(jù)的質(zhì)量直接影響下游分析的準(zhǔn)確性，而不同測(cè)序平臺(tái)的錯(cuò)誤率分布和偏差特性存在顯著差異阅羹。因此勺疼，在分析測(cè)序數(shù)據(jù)前，必須清晰了解以下幾點(diǎn)：
（1）測(cè)序平臺(tái)與錯(cuò)誤率分布
-?確定數(shù)據(jù)來(lái)源的測(cè)序平臺(tái)及其錯(cuò)誤率分布特征捏鱼。
- 判斷測(cè)序數(shù)據(jù)是否存在偏向性或局限性执庐，例如是否顯著受 GC 含量的影響。
GC含量是指DNA或RNA序列中鳥(niǎo)嘌呤（G）和胞嘧啶（C）這兩種核苷酸的比例导梆。GC含量通常用以下公式計(jì)算：
- 了解平臺(tái)特性的同時(shí)，結(jié)合官方文檔與自身數(shù)據(jù)分析進(jìn)行驗(yàn)證看尼。
（2） 錯(cuò)誤率對(duì)下游分析的潛在影響
- 變異檢測(cè)：測(cè)序深度是否充足递鹉？變異位點(diǎn)是否過(guò)度集中在 read 的末端？是否存在正反鏈偏向性藏斩？
- 基因組組裝：組裝過(guò)程中對(duì)錯(cuò)誤堿基的敏感性較高躏结，需要更嚴(yán)格的剪切策略，以避免計(jì)算資源和時(shí)間的額外消耗狰域。
（3）實(shí)際策略與建議
- 全面了解數(shù)據(jù)：分析數(shù)據(jù)質(zhì)量圖譜媳拴，評(píng)估錯(cuò)誤率分布和 GC 偏向性黄橘，必要時(shí)進(jìn)行質(zhì)量剪切和糾錯(cuò)。
- 視具體需求調(diào)整策略：如變異分析中評(píng)估深度與偏向性屈溉，或基因組組裝時(shí)優(yōu)化數(shù)據(jù)預(yù)處理流程塞关。
二、如何全面認(rèn)識(shí)原始測(cè)序數(shù)據(jù)（FASTQ 數(shù)據(jù)）
原始測(cè)序數(shù)據(jù)的分析是了解實(shí)驗(yàn)質(zhì)量和評(píng)估下游分析潛在問(wèn)題的第一步子巾。一般來(lái)說(shuō)帆赢，我們可以從以下幾個(gè)方面入手：
（1）Read 各個(gè)位置的堿基質(zhì)量值分布
通過(guò)評(píng)估每個(gè)位置上的質(zhì)量值，可以直觀反映測(cè)序錯(cuò)誤率的趨勢(shì)线梗。
（2）堿基的總體質(zhì)量值分布
分析所有 read 的質(zhì)量值概覽椰于，了解數(shù)據(jù)的整體可靠性。
（3）Read 各位置的堿基分布比例
目的是分析各堿基 (A/T/G/C) 的分離程度是否正常缠导，例如是否存在某些堿基過(guò)多或過(guò)少的現(xiàn)象廉羔。
（4）GC 含量分布
檢查測(cè)序數(shù)據(jù)是否存在 GC 偏向性或異常波動(dòng)溉痢，這可能會(huì)影響下游分析僻造。
（5）Read 各位置的 N 含量
N 表示不確定堿基，其分布能反映測(cè)序錯(cuò)誤的隨機(jī)性或局部錯(cuò)誤的集中性孩饼。
（6）Read 是否包含測(cè)序接頭序列
接頭序列的存在會(huì)干擾比對(duì)和組裝髓削，需檢測(cè)并剪切。
（7）Read 重復(fù)率
評(píng)估因 PCR 擴(kuò)增引入的重復(fù)情況镀娶，重復(fù)率過(guò)高可能表明樣本復(fù)雜度不足或?qū)嶒?yàn)問(wèn)題立膛。
通過(guò)全面分析以上幾個(gè)方面，可以對(duì)原始數(shù)據(jù)的基本質(zhì)量和潛在問(wèn)題有一個(gè)清晰的認(rèn)識(shí)梯码。
三宝泵、Read 各位置堿基質(zhì)量值分布的分析
在前文中已經(jīng)展示了如何使用 Python 代碼計(jì)算單條 read 的質(zhì)量值和測(cè)序錯(cuò)誤率。然而轩娶，當(dāng)面對(duì)成千上萬(wàn)條 read 時(shí)儿奶，逐條分析顯然不現(xiàn)實(shí)，此時(shí)更直觀的表達(dá)方式是數(shù)據(jù)可視化鳄抒。
FastQC是目前廣泛使用的工具之一闯捎，它能夠高效生成測(cè)序數(shù)據(jù)的質(zhì)量控制 (QC) 報(bào)告。其報(bào)告涵蓋了上述多個(gè)方面许溅，并提示數(shù)據(jù)可能存在的問(wèn)題瓤鼻。例如，F(xiàn)astQC 會(huì)通過(guò)綜合分析 FastQ 文件中所有 read 的數(shù)據(jù)贤重，為 read 各位置的堿基質(zhì)量值分布繪制圖表茬祷，如下圖所示：
注：在這個(gè)分布圖中，橫軸表示 read 的位置并蝗，縱軸表示質(zhì)量值祭犯。對(duì)于一個(gè)測(cè)序質(zhì)量較好的數(shù)據(jù)集耐朴，大多數(shù)堿基的質(zhì)量值應(yīng)在高分區(qū)域（通常為綠色）。FastQC 的優(yōu)勢(shì)在于快速生成報(bào)告并幫助識(shí)別潛在問(wèn)題盹憎，但其缺點(diǎn)是圖表美觀度較差啦~
除了上面read各位置的堿基質(zhì)量值分布之外筛峭，F(xiàn)astQC還會(huì)為我們計(jì)算其他幾個(gè)非常有價(jià)值的統(tǒng)計(jì)結(jié)果，包括：
1）堿基總體質(zhì)量值分布陪每，只要大部分都高于20影晓，那么就比較正常。
2）read各個(gè)位置上堿基比例分布
3）GC含量分布圖
4）N含量分布圖
5）接頭序列
（1）測(cè)序接頭與其在數(shù)據(jù)中的影響
在第 1 節(jié)中提到，測(cè)序文庫(kù)構(gòu)建是測(cè)序流程中的關(guān)鍵一步对湃，而**測(cè)序接頭（adapter）**正是在這個(gè)過(guò)程中被添加的崖叫。接頭的作用有以下兩點(diǎn)：
① 固定與區(qū)分樣本
接頭用于將 DNA 片段固定到測(cè)序平臺(tái)的流動(dòng)槽（flowcell）上。
在多樣本混合測(cè)序時(shí)拍柒，接頭上還包含特定的條形碼序列（barcode）心傀，以便在測(cè)序后區(qū)分不同樣本。
② 連接與擴(kuò)增
接頭為插入片段的末端提供連接點(diǎn)拆讯，便于后續(xù)擴(kuò)增和測(cè)序脂男。
（2）接頭序列的測(cè)到與其影響
當(dāng)測(cè)序 read 的長(zhǎng)度大于目標(biāo) DNA 片段（插入片段）的長(zhǎng)度時(shí)，read 的末尾可能會(huì)延伸到接頭序列种呐。
① WGS（全基因組測(cè)序）中的情況
在 WGS 測(cè)序中宰翅，這種現(xiàn)象較少發(fā)生。構(gòu)建 WGS 文庫(kù)時(shí)爽室，插入片段通常較長(zhǎng)（幾百 bp）汁讼，而測(cè)序 read 的長(zhǎng)度一般為 100–150 bp。因此阔墩，read 通常不會(huì)測(cè)到接頭序列嘿架。
② RNA 測(cè)序中的情況
在一些 RNA 測(cè)序（如 miRNA 測(cè)序）中，由于目標(biāo) RNA 片段本身較短（通常幾十 bp）戈擒，read 很容易測(cè)完整個(gè)片段并延伸到接頭序列眶明。這種現(xiàn)象在短片段 RNA 測(cè)序中尤為常見(jiàn)。
測(cè)到的接頭序列會(huì)干擾下游分析筐高，其影響類似于低質(zhì)量堿基。未切除的接頭序列可能導(dǎo)致：
比對(duì)錯(cuò)誤丑瞧，降低比對(duì)率和精確性柑土；
錯(cuò)誤的變異檢測(cè)或基因表達(dá)量估計(jì)。
（3）接頭序列的清理
在正式數(shù)據(jù)分析前绊汹，接頭序列與低質(zhì)量堿基一樣稽屏，需要通過(guò)預(yù)處理軟件清理。常用的清理工具包括：
① Trimmomatic：高效去除接頭序列與低質(zhì)量片段西乖；
② Cutadapt：能夠快速識(shí)別并剪切指定的接頭序列狐榔；
③ FASTP：支持自動(dòng)檢測(cè)接頭序列并進(jìn)行質(zhì)量控制。
清理后的數(shù)據(jù)可以顯著提高下游分析的準(zhǔn)確性获雕，避免因接頭干擾而消耗計(jì)算資源薄腻。
通過(guò)預(yù)處理確保測(cè)序數(shù)據(jù)的高質(zhì)量是后續(xù)研究的基礎(chǔ)，也是數(shù)據(jù)可靠性的保障届案。
四庵楷、FastQC該怎么用呢？
（1）通過(guò)網(wǎng)頁(yè)搜索或者直接到它的主頁(yè)上下載最新的版本。
（2）FastQC的運(yùn)行非常簡(jiǎn)單，直接在終端通過(guò)命令行是最有效直接的弄贿，例如：
$ /path_to_fastqc/FastQC/fastqc untreated.fq -o fastqc_out_dir/
注意：-o 參數(shù)用于指定FastQC報(bào)告的輸出目錄春锋，這個(gè)目錄需要事先創(chuàng)建好，如果不指定特定的目錄差凹，那么FastQC的結(jié)果會(huì)默認(rèn)輸出到文件untreated.fq的同一個(gè)目錄下看疙。它輸出結(jié)果只有兩個(gè)，一個(gè)html和一個(gè).zip壓縮包直奋。
$ tree fastqc_out_dir/
（3）我們可以直接通過(guò)瀏覽器打開(kāi)html能庆，就可以看到FastQC給出的所有結(jié)果，zip壓縮包解壓后脚线，從中我們也可以在對(duì)應(yīng)的目錄下找到所有的QC圖表和Summary數(shù)據(jù)搁胆。
除了上述用法之外，F(xiàn)astQC支持同時(shí)輸入多個(gè)fq文件（或者以通配符的形式輸入fq）邮绿，當(dāng)我們的fq文件比較多時(shí)渠旁，這種用法會(huì)比較方便，如：
$ /path_to_fastqc/FastQC/fastqc /path_to_fq/*.fq -o fastqc_out_dir/
這個(gè)鏈接是FastQC的一個(gè)結(jié)果模板：Online report
五船逮、切除測(cè)序接頭序列和read的低質(zhì)量序列
（1）接頭序列與低質(zhì)量堿基的切除工具及原理
目前市面上已有多種工具用于切除測(cè)序接頭序列和低質(zhì)量堿基顾腊，其中主流工具包括SOAPnuke、cutadapt挖胃、untrimmed等十余種杂靶，但在便捷性和功能性上，以下幾種工具尤為突出：
① Trimmomatic：功能全面酱鸭，支持接頭序列切除與低質(zhì)量堿基過(guò)濾吗垮。
② Sickle：專注于低質(zhì)量堿基過(guò)濾，支持雙端（PE）和單端（SE）數(shù)據(jù)凹髓。
③ Seqtk：快速處理低質(zhì)量堿基烁登，但僅支持單端（SE）數(shù)據(jù)。
（2）Trimmomatic 的優(yōu)勢(shì)
Trimmomatic不僅可以切除 Illumina 平臺(tái)默認(rèn)的接頭序列蔚舀，還支持切除用戶自定義的接頭序列饵沧，同時(shí)過(guò)濾 read 末尾的低質(zhì)量堿基。相比之下赌躺，Sickle和Seqtk僅用于去除低質(zhì)量片段狼牺，不支持接頭切除。
① 工具工作原理
以Trimmomatic為例寿谴，其工作原理基于滑動(dòng)窗口法：
- 在指定窗口長(zhǎng)度內(nèi)锁右，計(jì)算窗口內(nèi)堿基的平均質(zhì)量值。
- 如果窗口的平均質(zhì)量值低于閾值，從該窗口開(kāi)始咏瑟，將 read 的尾部全部切除拂到，而不是挖掉局部低質(zhì)量片段。
- 這種方式保證了 read 的連續(xù)性码泞，避免人為改變實(shí)際的基因組序列兄旬。
特別注意：切除后的序列可能會(huì)顯著縮短，因此在實(shí)際操作中需要設(shè)定最低 read 長(zhǎng)度以避免無(wú)效數(shù)據(jù)進(jìn)入下游分析余寥。
② 工具選擇建議
- 如果下機(jī)數(shù)據(jù)中包含接頭序列领铐，推薦使用Trimmomatic，它功能全面宋舷，且對(duì)多種場(chǎng)景適用绪撵。
- 如果不包含接頭序列，且對(duì)計(jì)算資源要求較高祝蝠，推薦使用Sickle（支持 PE 和 SE）或Seqtk（僅支持 SE）音诈，兩者處理速度更快且資源占用較低。
(3) 使用 Trimmomatic 過(guò)濾數(shù)據(jù)
安裝與準(zhǔn)備
1. 獲取 Trimmomatic
- 最新版本可從其官網(wǎng)下載绎狭。
- 如需查看源代碼细溅，可在GitHub上找到項(xiàng)目頁(yè)面。
2. 安裝與運(yùn)行
- 下載打包的二進(jìn)制文件（如 Trimmomatic-0.36.zip）儡嘶，解壓后可看到trimmomatic-*.jar文件喇聊。
- 使用 Java 運(yùn)行程序：java -jar trimmomatic-0.36.jar
3. 構(gòu)建過(guò)濾流程
使用 Trimmomatic 時(shí)，通常通過(guò)命令行指定輸入文件蹦狂、接頭序列文件和過(guò)濾參數(shù)誓篱。以下是一個(gè)基本流程示例：
ILLUMINACLIP：指定接頭序列文件及參數(shù)（最大錯(cuò)配數(shù)、精確性評(píng)分等）鸥咖。
SLIDINGWINDOW：設(shè)置滑動(dòng)窗口長(zhǎng)度和質(zhì)量閾值燕鸽。
MINLEN：過(guò)濾掉低于指定長(zhǎng)度的序列。
4.輸出文件解釋
output_R1_paired.fastq.gz和output_R2_paired.fastq.gz：過(guò)濾后仍成對(duì)的 reads啼辣。
output_R1_unpaired.fastq.gz和output_R2_unpaired.fastq.gz：未成對(duì)的 reads，通常在后續(xù)分析中被丟棄御滩。
接頭序列的選擇與設(shè)置
（1）接頭序列的選擇
① 接頭序列
- HiSeq和MiSeq系列通常使用TruSeq3接頭序列鸥拧。
- TruSeq2接頭序列主要用于較老的 GA2 測(cè)序儀，目前已較少使用削解。
- 這些默認(rèn)接頭序列信息可以在 Illumina 官網(wǎng)查詢富弦。
② 自定義接頭序列
- 如果所用測(cè)序平臺(tái)非 Illumina，或者需要切除自定義接頭序列氛驮，可以自行創(chuàng)建符合格式的接頭序列文件腕柜，并作為參數(shù)傳入 Trimmomatic。
- 具體命名和格式要求可參考 Trimmomatic 文檔（The Adapter Fasta）。
（2）PE 與 SE 模式選擇
① 接頭序列的使用需根據(jù)具體測(cè)序類型選擇：
- **PE 模式（Pair End）：**雙端測(cè)序盏缤，使用TruSeq3-PE.fa文件砰蠢。
- **SE 模式（Single End）：**單端測(cè)序，使用TruSeq3-SE.fa文件唉铜。
Trimmomatic 的運(yùn)行模式（兩種）
① PE 模式：適用于雙端測(cè)序數(shù)據(jù)台舱。
PE 模式命令：java -jar trimmomatic-0.36.jar PE [參數(shù)] 輸入文件 輸出文件 過(guò)濾參數(shù)
示例命令
PE 模式（HiSeq PE 測(cè)序）
② SE 模式：適用于單端測(cè)序數(shù)據(jù)。
SE 模式命令：java -jar trimmomatic-0.36.jar SE [參數(shù)] 輸入文件 輸出文件 過(guò)濾參數(shù)
SE 模式（HiSeq SE 測(cè)序）
參數(shù)說(shuō)明
常用參數(shù)
- phred33 或 -phred64?指定堿基質(zhì)量值體系潭流。
?- Phred33是當(dāng)前測(cè)序數(shù)據(jù)的主流格式竞惋，必須顯式指定此參數(shù)。
?- 若未指定灰嫉，默認(rèn)使用Phred64拆宛，可能導(dǎo)致錯(cuò)誤處理。
- trimlog?指定修剪日志文件讼撒，記錄每個(gè) read 的修剪信息浑厚。
- ILLUMINACLIP
指定接頭序列文件和參數(shù)，格式為：ILLUMINACLIP:<文件路徑>:<最大錯(cuò)配數(shù)>:<精確性閾值>:<剪切長(zhǎng)度>
- SLIDINGWINDOW
滑動(dòng)窗口參數(shù)椿肩，格式為：SLIDINGWINDOW:<窗口大小>:<質(zhì)量閾值>
- LEADING 和 TRAILING?分別指定去除 read 起始或末端低質(zhì)量堿基的質(zhì)量閾值瞻颂。
- MINLEN?設(shè)置修剪后 read 的最低長(zhǎng)度，低于此值的 read 將被丟棄郑象。
PE 與 SE 模式的差異
PE 模式
- 需要輸出過(guò)濾后和未過(guò)濾的雙端 read 數(shù)據(jù)：
- out.read_1.fq.gz和out.read_2.fq.gz：過(guò)濾后的配對(duì) reads贡这。
- out.trim.read_1.fq.gz和out.trim.read_2.fq.gz：未通過(guò)過(guò)濾的 reads。
SE 模式
- 不需要指定未通過(guò)過(guò)濾的 reads 文件厂榛。過(guò)濾后直接生成一個(gè)輸出文件盖矫。
- 例如，在 SE 模式下：raw_data/untreated.fq --> out.untreated.fq.gz
- 通過(guò)靈活設(shè)置接頭序列和運(yùn)行模式击奶，Trimmomatic可以高效完成測(cè)序數(shù)據(jù)的預(yù)處理辈双，為下游分析提供高質(zhì)量的輸入數(shù)據(jù)。
以下是對(duì)Trimmomatic參數(shù)的整理：
ILLUMINACLIP 參數(shù)解釋
- 接頭序列文件：指定用于匹配的接頭序列文件（如TruSeq3-PE.fa或自定義文件）柜砾。
- 最大錯(cuò)配數(shù)：允許接頭序列匹配時(shí)的最大錯(cuò)配數(shù)湃望，值越小匹配越嚴(yán)格。
- 雙端 reads 匹配度：兩條 reads 同時(shí)與接頭序列匹配時(shí)痰驱，總體匹配度的最低閾值（推薦值為 30%）证芭。
- 單端匹配度：?jiǎn)螚l read 與接頭序列部分匹配時(shí)的最低匹配度（推薦值為 10%）。
**注：**實(shí)際測(cè)序數(shù)據(jù)可能只覆蓋接頭序列的一部分担映，因此不要求完全匹配废士。上述推薦值（30% 和 10%）是常用參數(shù)。
以上就是數(shù)據(jù)質(zhì)控的內(nèi)容結(jié)束蝇完。
生物信息學(xué)領(lǐng)域非常廣泛官硝，難以一次說(shuō)盡矗蕊。我們下次繼續(xù)更新，一起深入學(xué)習(xí)生物信息學(xué)的內(nèi)容氢架！
喜歡的寶子們點(diǎn)個(gè)贊吧~碼字不易傻咖，且行且珍惜~
