前言：
????一般測(cè)序下機(jī)數(shù)據(jù)會(huì)存在含N比例過大、測(cè)序質(zhì)量較低的堿基數(shù)占比過高、含有duplication、序列污染等低質(zhì)量reads,這些不合格的reads會(huì)影響后續(xù)的分析，所以，我們拿到測(cè)序數(shù)據(jù)首先要了解測(cè)序數(shù)據(jù)的質(zhì)量情況，具體內(nèi)容包括含N比例、GC含量、duplication情況、序列長(zhǎng)度分布情況、堿基平衡情況等。
????今天，我們將一起通過數(shù)據(jù)格式和質(zhì)量體系、數(shù)據(jù)質(zhì)控步驟、Fastqc結(jié)果解讀及異常處理三大模塊進(jìn)行學(xué)習(xí)。
第一部分?數(shù)據(jù)格式和質(zhì)量體系
????Illumina測(cè)序的下機(jī)數(shù)據(jù)一般為fastq格式，至于fastq格式的說明我已經(jīng)在上期《測(cè)序技術(shù)原理及常用數(shù)據(jù)格式簡(jiǎn)介》中有詳細(xì)描述，在進(jìn)行數(shù)據(jù)質(zhì)控前，我們需要知道數(shù)據(jù)中第四行質(zhì)量字符和序列質(zhì)量值Q值的關(guān)系以及Q值與堿基測(cè)序錯(cuò)誤率的關(guān)系。
????Fastq數(shù)據(jù)中的質(zhì)量字符并不是和質(zhì)量值Q值直接對(duì)應(yīng)起來的，fastq數(shù)據(jù)格式中的質(zhì)量字符是ASCII值，在Phred+64體系中，ASCII值-64的結(jié)果就是Q值，在Phred+33體系中，ASCII值-33的結(jié)果就是Q值。在Phred+33體系中，Q = -10log10(P), 堿基質(zhì)量值與誤率的對(duì)應(yīng)關(guān)系表如下所示：
即，Q10準(zhǔn)確率為90%，Q20準(zhǔn)確率為99%，Q30準(zhǔn)確率為99.9%，Q40準(zhǔn)確率為99.99%，Q50準(zhǔn)確率為99.999%。
第二部分 數(shù)據(jù)質(zhì)控
????數(shù)據(jù)質(zhì)控現(xiàn)在用得最多的是fastqc,我們今天就以它為工具學(xué)習(xí)如何了解測(cè)序數(shù)據(jù)質(zhì)量。
Fastqc下載安裝
wget -c http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.5.zip
unzip fastqc_v0.11.5.zip
cd FastQC
chmod +x fastqc
Fastqc評(píng)估測(cè)序數(shù)據(jù)質(zhì)量
Usage:
fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam] [-c contaminant file][-t]
--(no)extract輸出的結(jié)果不接壓，若無此選項(xiàng)，輸出的結(jié)果為.zip壓縮文件。
-f fastq|bam|sam指定輸入文件格式，若無此項(xiàng)，則會(huì)自動(dòng)檢測(cè)。
-c contaminant?file指定一個(gè)contaminant文件，文件格式為”Name\tSequence”，fastqc會(huì)把overrepreseted sequence往這個(gè)contaminant文件搜索。
-t線程數(shù)
例子：
fastqc *fq.gz –t 4 ?#目錄下所有fq.gz文件進(jìn)行質(zhì)控，線程數(shù)一般與文件數(shù)一致。
第三部分?fastqc結(jié)果解讀及異常處理
????橫軸代表堿基在序列中的位置，縱軸代表Q值，由前面堿基質(zhì)量值與錯(cuò)誤率的關(guān)系可知，若某個(gè)位置對(duì)應(yīng)的Q值為30，則該處堿基測(cè)序準(zhǔn)確率為99.9%。
????如Figure1所示，在箱線圖中，紅色表示中位數(shù)，黃色是25%-75%區(qū)間，觸須是10%-90%區(qū)間，藍(lán)線是平均數(shù)。若任一位置的下四分位數(shù)低于10或中位數(shù)低于25，報(bào)"WARN"；若任一位置的下四分位數(shù)低于5或中位數(shù)低于20，報(bào)"FAIL"。
????當(dāng)出現(xiàn)任一位置的下四分位數(shù)低于10或中位數(shù)低于25或任一位置的下四分位數(shù)低于5或中位數(shù)低于20時(shí)，表示測(cè)序數(shù)據(jù)存在質(zhì)量不合格的情況，這時(shí)我們可以繼續(xù)觀察Sequence Contentacross圖、GC Contentacross all base圖、N Content across all bases圖、Sequences?Duplication level stastics圖這幾個(gè)圖進(jìn)一步判斷測(cè)序數(shù)據(jù)的不合格之處具體在哪。
????觀察Sequence Contentacross圖和GC Contentacross all base圖的GC含量的線是否平行于X軸，若不平行，則該位置往往有overrepresentedsequence的污染，可能原因建庫(kù)過程的誤差、測(cè)序的系統(tǒng)誤差或者文庫(kù)本身特點(diǎn)。
????由N Content across all bases圖可知reads中含N堿基的情況，理想狀況下是含N量越少越好，在微生物多樣性分析中一般是去除含N堿基比例>5%的序列。
????去除含N比例過高的序列，可以用NGS QC tookit。該工具可在http://www.nipgr.res.in/ngsqctoolkit.html處下載，解壓縮在QC/RIMINGREADS文件夾中用AmbiguityFiltering.pl腳本去除含N比例過高的序列，具體使用方法可以參考軟件壓縮包內(nèi)的manual文件。
????觀察Sequences?Duplication level stastics圖，橫坐標(biāo)是duplication的次數(shù)，縱坐標(biāo)是duplicated reads的數(shù)目，若duplication的程度偏高，則可能存在PCR duplication。去除duplication可以通過Samtools、Picard或Iontorrent，其中Samtools只看5’端的起始位置不考慮reads突變；Picard不僅考慮起始位點(diǎn)也會(huì)考慮突變情況和質(zhì)量值，即reads完全一樣的才會(huì)被當(dāng)成duplication被去除；Iontorrent則是看5’端的起始位置和3’端adaptor的比對(duì)情況，不考慮reads突變。
????去除PCR duplication的具體操作可參考：仔細(xì)探究samtools的rmdup是如何行使去除PCR重復(fù)reads功能的。
????最后，在進(jìn)行去低質(zhì)量reads和接頭等預(yù)處理步驟后，再次進(jìn)行fastqc質(zhì)控，然后用Multic QC即可把多個(gè)樣品的質(zhì)控結(jié)果匯總在一起，在報(bào)告中圖片是交互式的，鼠標(biāo)停留可顯示樣品名。如下圖所示。
????當(dāng)然，要是不想經(jīng)歷fastqc、NGS QC tookit、cutadapter、trimomatic這些工具這么麻煩，想省事點(diǎn)，直接一步到位完成所有測(cè)序數(shù)據(jù)預(yù)處理步驟，有沒有這樣的神器呢？有的，fastp可以做到。常規(guī)測(cè)序預(yù)處理流程一般包括Fastqc等軟件進(jìn)行質(zhì)控、cutadapter去除接頭、Trimmonatic進(jìn)行過濾低質(zhì)量reads、切除長(zhǎng)度不合格的序列等步驟，涉及軟件較多，比較繁瑣，直接用fastp則可以一步到位完成。詳情請(qǐng)點(diǎn)擊https://github.com/OpenGene/fastp自行了解。
