20160410 測序分析——使用 FastQC 做質控

01-26

前言

Hello 大家好！我是孟浩巍，一名生物信息方向的搬磚工！

之前我們已經介紹過了二代測序Illumina平台的測序原理，包括一些常用的專業術語，比如reads，flowcell，tail等等，如果大家有什麼疑問請看之前的文章：

Illumina測序原理介紹

隨後我們又介紹了，測序數據的儲存格式FASTQ格式。今天要講的內容會涉及到一些裡面的知識。如果你忘記啦，請看看之前的文章：

FASTA與FASTQ格式介紹

Ok！那麼今天我們要解決的問題是在從測序公司拿到原始數據以後，我們應該怎麼評價這次的測序質量。是不是要做相應的一些後續處理。我們今天要使用的就是一個強大的工具——FastQC

FastQC的基本介紹

FastQC是一款基於Java的軟體，一般都是在linux環境下使用命令行運行，它可以快速多線程地對測序數據進行質量評估（Quality Control），其官網地址為：Babraham Bioinformatics

FastQC的下載和安裝，和一般的Java軟體沒有什麼區別，我們在這裡就不做介紹了，在成功安裝好以後，我們就在命令行模式下，輸入fastqc就可以調用這個程序，界面如下：

這時候我們可以選擇 --help選項看一下幫助文檔：

# 基本格式# fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam] [-c contaminant file] seqfile1 .. seqfileN# 主要是包括前面的各種選項和最後面的可以加入N個文件# -o --outdir FastQC生成的報告文件的儲存路徑，生成的報告的文件名是根據輸入來定的# --extract 生成的報告默認會打包成1個壓縮文件，使用這個參數是讓程序不打包# -t --threads 選擇程序運行的線程數，每個線程會佔用250MB內存，越多越快咯# -c --contaminants 污染物選項，輸入的是一個文件，格式是Name [Tab] Sequence，裡面是可能的污染序列，如果有這個選項，FastQC會在計算時候評估污染的情況，並在統計的時候進行分析，一般用不到# -a --adapters 也是輸入一個文件，文件的格式Name [Tab] Sequence，儲存的是測序的adpater序列信息，如果不輸入，目前版本的FastQC就按照通用引物來評估序列時候有adapter的殘留# -q --quiet 安靜運行模式，一般不選這個選項的時候，程序會實時報告運行的狀況。

以我平時用的一個真實的例子：

fastqc -o ./tmp.result/fastQC/ -t 6 ./tmp.data/fastq/H1EScell-dnase-2014-GSE56869_20151208_SRR1248176_1.fq

使用的數據是2014年Dnase Hi-C的測序數據，數據下載地址：

http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM1370433

運行一段時間以後，就會出現報告：

H1EScell-dnase-2014-GSE56869_20151208_SRR1248176_1.fq_fastqc.htmlH1EScell-dnase-2014-GSE56869_20151208_SRR1248176_1.fq_fastqc.zip

使用瀏覽器打開後綴是html的文件，就是圖表化的fastqc報告：

FastQC的報告介紹

總結信息

上圖中Summary的部分就是整個報告的目錄，整個報告分成若干個部分。合格會有個綠色的對勾，警告是個「!」，不合格是個紅色的叉子。

基本信息

# Encoding指測序平台的版本和相應的編碼版本號，這個在計算Phred反推error P的時候有用，如果不明白可以參考之前的文章。# Total Sequences記錄了輸入文本的reads的數量# Sequence length 是測序的長度# %GC 是我們需要重點關注的一個指標，這個值表示的是整體序列中的GC含量，這個數值一般是物種特意的，比如人類細胞就是42%左右。

序列測序質量統計

# 此圖中的橫軸是測序序列第1個鹼基到第101個鹼基# 縱軸是質量得分，Q = -10*log10（error P）即20表示1%的錯誤率，30表示0.1%# 圖中每1個boxplot，都是該位置的所有序列的測序質量的一個統計，上面的bar是90%分位數，下面的bar是10%分位數，箱子的中間的橫線是50%分位數，箱子的上邊是75%分位數，下邊是25%分位數# 圖中藍色的細線是各個位置的平均值的連線# 一般要求此圖中，所有位置的10%分位數大於20,也就是我們常說的Q20過濾# 所以上面的這個測序結果，需要把後面的87bp以後的序列切除，從而保證後續分析的正確性# Warning 報警如果任何鹼基質量低於10,或者是任何中位數低於25# Failure 報錯如果任何鹼基質量低於5,或者是任何中位數低於20

每個tail測序的情況

# 橫軸和之前一樣，代表101個鹼基的每個不同位置# 縱軸是tail的Index編號# 這個圖主要是為了防止，在測序過程中，某些tail受到不可控因素的影響而出現測序質量偏低# 藍色代表測序質量很高，暖色代表測序質量不高，如果某些tail出現暖色，可以在後續分析中把該tail測序的結果全部都去除

每條序列的測序質量統計

# 假如我測的1條序列長度為101bp，那麼這101個位置每個位置Q之的平均值就是這條reads的質量值# 該圖橫軸是0-40，表示Q值# 縱軸是每個值對應的reads數目# 我們的數據中，測序結果主要集中在高分中，證明測序質量良好！

GC 含量統計

# 橫軸是1 - 101 bp；縱軸是百分比# 圖中四條線代表A T C G在每個位置平均含量# 理論上來說，A和T應該相等，G和C應該相等，但是一般測序的時候，剛開始測序儀狀態不穩定，很可能出現上圖的情況。像這種情況，即使測序的得分很高，也需要cut開始部分的序列信息，一般像我碰到這種情況，會cut前面5bp

序列平均GC含量分布圖

# 橫軸是0 - 100%；縱軸是每條序列GC含量對應的數量# 藍色的線是程序根據經驗分布給出的理論值，紅色是真實值，兩個應該比較接近才比較好# 當紅色的線出現雙峰，基本肯定是混入了其他物種的DNA序列# 這張圖中的信息良好

序列測序長度統計

# 每次測序儀測出來的長度在理論上應該是完全相等的，但是總會有一些偏差# 比如此圖中，101bp是主要的，但是還是有少量的100和102bp的長度，不過數量比較少，不影響後續分析# 當測序的長度不同時，如果很嚴重，則表明測序儀在此次測序過程中產生的數據不可信

序列Adapter

# 此圖衡量的是序列中兩端adapter的情況# 如果在當時fastqc分析的時候-a選項沒有內容，則默認使用圖例中的四種通用adapter序列進行統計# 本例中adapter都已經去除，如果有adapter序列沒有去除乾淨的情況，在後續分析的時候需要先使用cutadapt軟體進行去接頭，這個軟體以後我會介紹

重複短序列

# 這個圖統計的是，在序列中某些特徵的短序列重複出現的次數# 我們可以看到1-8bp的時候圖例中的幾種短序列都出現了非常多的次數，一般來說，出現這種情況，要麼是adapter沒有去除乾淨，而又沒有使用-a參數；要麼就是序列本身可能重複度比較高，如建庫PCR的時候出現了bias# 對於這種情況，我的辦法是可以cut掉前面的一些長度，可以試著cut 5~8bp

本期總結與下期預告

在本片文章中，我們已經介紹了測序質量評估最常用的FastQC軟體，並詳細解讀了報告內容。希望能夠對大家有所幫助。可是我們還是留有一些問題，比如我一直說cut序列，去adapter，卻沒有給出工具的用法。所以下期我們的主題主要是，圍繞著FastQC中出現的問題使用不同的fastx_trimmer，cutadapt等工具進行修正。

希望各位多多點贊支持！

孟浩巍

——————————————————————————

另外歡迎各位參加我們的知乎Live：

1. 知乎Live：如何快速入門生物信息學（涉及內容：測序原理，生物信息學發展歷史，軟體的安裝與調試，入門路線圖，介紹了RNA-Seq的分析流程並給出實踐代碼）；

如何快速入門生物信息學

2. 知乎Live: 生信進階第1課-重複Nature文章 (涉及內容：肺癌相關研究現狀，RNA-Seq單細胞測序，RNA-Seq的建庫方法，RNA-Seq的分析流程細節，相關生信圖的繪製）；

生信進階第1課-重複Nature文章

3. 知乎Live：生信進階第2課-基因組序列

(涉及內容：介紹基因組的序列結構，hg19與hg38的區別，ENCODE計劃，常用的表觀組學實驗原理ChIP-Seq，Hi-C等，ChIP-Seq的標準處理流程，繪圖原理)

生信進階第2課-基因組序列

4. 知乎Live：不用編程怎麼做生物信息學

(涉及內容：介紹生物信息學入門的幾個層次，從命令行到圖形界面再到命令行，繪製生物進化樹，圖形界面分析平台，使用圖形界面處理RNA-Seq數據，使用圖形界面分析ChIP-Seq數據，UCSC genome browser，WashU genome browser)

不用編程怎麼做生物信息學