linux质控命令,质控软件 fastp
1引言
各位做生信的小伙伴都知道,對于下機的FASTQ數據需要進行質控和預處理,以保證下游分析輸入的數據都是干凈可靠的。通常我們都是使用FASTQC等軟件進行質控,使用cutadapt軟件去除接頭,使用Trimmomatic等軟件進行剪裁,然后使用一些自已開發的腳本進行過濾。這一過程可能需要使用多個軟件,相當繁瑣,而且速度較慢,這些軟件大多又不支持多線程,遇到較大的FASTQ文件,處理起來可真是讓人等得心急如焚。
所以今天給大家介紹一款新的軟件,名字叫fastp,它可以僅僅掃描FASTQ文件一次,就完成比FASTQC
cutadapt + Trimmomatic 這三個軟件加起來的功能還多很多的功能,而且速度上比僅僅使用Trimmomatic一個軟件還要快3倍左右,因為它使用C++開發,處處使用了高效算法,而且完美支持多線程!正因為其強大的功能和飛快的速度,其github項目在第一個版本發布至今才一個月左右,就已經收到了50多個star,這在生信軟件小領域里面算是神速了!而且該軟件更新頻繁,至今已發布了14版,每一版都有新功能加入!該項目的github地址請戳:https://github.com/OpenGene/fastp
2功能特點
先來看一看該軟件的功能列表:
l 對數據自動進行全方位質控,生成人性化的報告。
l 過濾功能(低質量,太短,太多N……)。
l 對每一個序列的頭部或尾部,計算滑動窗內的質量均值,并將均值較低的子序列進行切除(類似Trimmomatic的做法,但是快非常多)。
l 全局剪裁 (在頭/尾部,不影響去重),對于Illumina下機數據往往最后一到兩個cycle需要這樣處理。
l 去除接頭污染。厲害的是,你不用輸入接頭序列,因為算法會自動識別接頭序列并進行剪裁。
l 對于雙端測序(PE)的數據,軟件會自動查找每一對read的重疊區域,并對該重疊區域中不匹配的堿基對進行校正。
l 去除尾部的polyG。對于Illumina NextSeq/NovaSeq的測序數據,因為是兩色法發光,polyG是常有的事,所以該特性對該兩類測序平臺默認打開。
l 對于PE數據中的overlap區間中不一致的堿基對,依據質量值進行校正
l 可以對帶分子標簽(UMI)的數據進行預處理,不管UMI在插入片段還是在index上,都可以輕松處理。
l 可以將輸出進行分拆,而且支持兩種模式,分別是指定分拆的個數,或者分拆后每個文件的行數。
以上功能大多都不需要輸入太多的參數,一些功能默認已經開啟,但是可以用參數關閉。fastp完美支持gzip的輸入和輸出,同時支持SE和PE數據,而且不但支持像Illumina平臺的short read數據,也在一定程度上支持了PacBio/Nanopore的long reads數據。
fastp軟件會生成HTML格式的報告,而且該報告中沒有任何一張靜態圖片,所有的圖表都是使用JavaScript動態繪制,非常具有交互性。想要看一下樣板報告的,可以去以下鏈接:http://opengene.org/fastp/fastp.html
而且軟件的開發者還充分考慮到了各種自動化分析的需求,不但生成了人可讀的HTML報告,還生成了程序可讀性非常強的JSON結果,該JSON報告中的數據包含了HTML報告100%的信息,而且該JSON文件的格式還是特殊定制的,不但程序讀得爽,你用任何一款文本編輯器打開,一眼過去也會看得明明白白。想要看一下JSON結果長什么樣的,可以去以下鏈接:http://opengene.org/fastp/fastp.json
3輕松上手
看起來這個軟件功能非常多,那使用起來是不是非常復雜呢?非也!該軟件的使用非常簡單,默認情況下只需要指定輸入和輸出文件,就可以很好地工作。例如我們想對in.fq文件進行過濾和質控,并輸出clean data為out.fq,那么就可以使用以下的命令:
fastp
-i in.fq -o out.fq
即使用小寫的i和小寫的o分別指定input和output文件,就大功告成啦。軟件執行完成之后,會生成out.fq,還會生成兩個文件fastp.html和fastp.json。其中fastp.html是可視化的QC質控報告以及各類過濾統計,而fastp.json是JSON版本的報告,主要用于下游程序來解讀質控和過濾的結果。而且fastp會同時統計過濾前(raw data)和過濾后(clean data)的質量信息,以方便你分析過濾前后數據質量發生了什么變化,夠貼心吧?
以上例子是對單端測序數據(single-end,SE)進行的,那對于雙端測序數據(paired-end,PE)是不是也可以呢?答案自然是肯定的,而且命令也基本上差不多,示例如下:
fastp
-i in.R1.fq -o out.R1.fq -I in.R2.fq -O out.R2.fq
可以看到,-i和-o還是用來指定read1的輸入了輸出,而大寫的-I和-O(注意是喔,而不是零)則是用于指定read2的輸入和輸出,其他都保持不變。而且fastp軟件最初的研發就是為了更好地處理PE數據,所以對于PE數據開發了更多的算法,比如基于overlap分析進行堿基校正等功能,就是只有PE數據獨享的。
fastp對于輸入和輸出都支持gzip壓縮,使用方法也很簡單,只要文件名的末尾帶有.gz,就會被認為是gzip壓縮文件,會啟用gzip對輸入輸出進行壓縮和解壓處理,例如以上PE的例子,如果是壓縮的,就可以是以下命令:
fastp
-i in.R1.fq.gz -I in.R2.fq.gz -o out.R1.fq.gz -O out.R2.fq.gz
4安裝fastp
安裝fastp十分簡單,如果你使用的是Linux系統,可以直接使用網站上提供的預編譯好的版本,下載地址是 http://opengene.org/fastp/fastp,可以使用wget等命令進行下載,下載了之后記得使用 chmod a+x ./fastp 增加該文件的可執行權限,然后就可以使用了。
也可以從源代碼進行編譯,需要使用git工具或者直接在網頁上下載release的源代碼,以git下載最新的代碼為例:
cd
fastp
make
sudo
make install
5功能介紹
接下來我們簡單介紹一下fastp的一些功能,受公眾號篇幅影響,每一個功能我們都只是簡單地帶過,如果想要看詳細的介紹,可以上github官網上查看。
5.1過濾
fastp可以對低質量序列,較多N的序列,該功能默認是啟用的,但可以使用-Q參數關閉。使用-q參數來指定合格的phred質量值,比如-q 15表示質量值大于等于Q15的即為合格,然后使用-u參數來指定最多可以有多少百分比的質量不合格堿基。比如-q 15 -u 40表示一個read最多只能有40%的堿基的質量值低于Q15,否則會被扔掉。使用-n可以限定一個read中最多能有多少個N。
fastp還默認啟用了read長度過濾,但可以使用-L參數關閉。使用-l參數指定最低要求一個read有多長,比如-l 30表示低于30個堿基的read會被扔掉。這個功能可以用于實現常用的discard模式,以保證所有輸出的序列都一樣長。
在fastp的HTML報告中,最頭上的Summary表格很清楚地顯示了過濾的統計信息,如下圖所示:
image
[圖片上傳失敗...(image-1b7d16-1576761117648)]
5.2接頭處理
接頭(adapter)污染的處理是FASTQ文件預處理中很重要的一步。fastp默認啟用了接頭處理,但是可以使用-A命令來關掉。fastp可以自動化地查找接頭序列并進行剪裁,也就是說你可以不輸入任何的接頭序列,fastp全自動搞定了!對于SE數據,你還是可以-a參數來輸入你的接頭,而對于PE數據則完全沒有必要,fastp基于PE數據的overlap分析可以更準確地查找接頭,去得更干凈,而且對于一些接頭本身就有堿基不匹配情況處理得更好。fastp對于接頭去除會有一個匯總的報告,如下圖所示:
image
[圖片上傳失敗...(image-1efdeb-1576761117648)]
5.3滑窗質量剪裁
很多時候,一個read的低質量序列都是集中在read的末端,也有少部分是在read的開頭。fastp支持像Trimmomatic那樣對滑動窗口中的堿基計算平均質量值,然后將不符合的滑窗直接剪裁掉。使用-5參數開啟在5’端,也就是read的開頭的剪裁,使用-3參數開啟在3’端,也就是read的末尾的剪裁。使用-W參數指定滑動窗大小,默認是4,使用-M參數指定要求的平均質量值,默認是20,也就是Q20。
5.4 PE數據的堿基校正
fastp支持對PE數據的每一對read進行分析,查找它們的overlap區間,然后對于overlap區間中不一致的堿基,如果發現其中一個質量非常高,而另一個非常低,則可以將非常低質量的堿基改為相應的非常高質量值的堿基值,如下圖所示:
image
[圖片上傳失敗...(image-bbf731-1576761117648)]
上圖中所示的標紅的T堿基是低質量序列,和高質量的A不匹配,它會被校正為A。該校正功能默認沒有開啟使用-c參數可以啟用,對于一些對噪聲容忍度低的應用,比如液體活檢,建議開啟。
5.5全局剪裁
fastp可以對所有read在頭部和尾部進行統一剪裁,該功能在去除一些測序質量不好的cycle比較有用,比如151*2的PE測序中,最后一個cycle通常質量是非常低的,需要剪裁掉。使用-f和-t分別指定read1的頭部和尾部的剪裁,使用-F和-T分別指定read2的頭部和尾部的剪裁。
5.6 polyG剪裁
對于兩色發光法的Illumina設備(NextSeq /
NovaSeq),因為在沒有光信號情況下base calling的結果會返回G,所以在序列的尾端可能會出現較多的polyG,需要被去除。fastp會自動化地識別NextSeq / NovaSeq的數據,然后進行polyG識別和剪裁。如果你想強制開啟該功能,可以指定-g參數,如果想強制關閉該功能,則可以指定-G參數。
5.7分子標簽UMI處理
UMI在處理ctDNA類似的超低頻突變檢測應用中是十分有用的,為了更好地對帶UMI的FASTQ文件進行預處理,fastp也很好地支持了UMI預處理功能。該功能默認沒有啟用,需要使用-U參數開啟,另外需要使用--umi_loc來指定UMI所在的位置,它可以是(index1、 index2、 read1、 read2、 per_index、 per_read )中的一種,分別表示UMI是在index位置上,還是在插入片段中。如果指定了是在插入序列中,還需要使用 --umi_len 參數來指定UMI所占的堿基長度。
5.8輸出文件切分
很多時候我們需要對輸出的FASTQ進行切分,分成大小均勻的多個文件,這樣可以使用比對軟件并行地比對,提高并行處理的速度。fastp軟件也提供了相應的功能,并且支持了兩種模式,分別是使用參數-s指定切分后文件的個數,或者使用-S參數指定每個切分后文件的行數。
6質控報告解讀
接下來,我們再看一下如何理解fastp生成的質控報告。fastp的報告在單一文件中同時包含了過濾前和過濾后的統計結果,如果是PE數據,則同時包含了read1和read2的統計結果。之前我們已經說過了,fastp會生成HTML的報告和JSON格式的報告。HTML報告的默認文件名是fastp.html,但是可以通過-h參數修改,JSON報告的默認文件名是fastp.json,但是可以通過-j參數修改。而且fastp報告還有一個標題,默認是fastp report,這個也可以通過-R參數修改為你想要的標題。JSON格式的報告是優化過的,人機皆可讀,適合進階的用戶使用程序解析,而這里我們重點關注HTML格式的報告。
6.1質量分布曲線圖
我們第一關注的當然是質量,所以fastp提供了質量分布曲線,即每一個cycle的平均質量值,而且fastp同時提供了A/T/C/G四種不同堿基的平均質量,以及總的平均質量,如下圖所示:
image
[圖片上傳失敗...(image-9b70f0-1576761117648)]
從上圖我們可以看到,一共有5條曲線,分別是A/T/C/G和mean。而且HTML報告中的每一個項目和分項目都是可以點擊進行隱藏和展開的。
6.2堿基含量分布曲線
和質量分布曲線類似,堿基含量分布曲線也是按照每一個cycle來的,顯示了每一個位置的堿基含量。如下圖所示:
image
[圖片上傳失敗...(image-6d6731-1576761117648)]
image
[圖片上傳失敗...(image-38c50f-1576761117648)]
從圖中可以看到,fastp同時顯示了A/T/C/G/N/GC
的每一個位置的比例和總的比例。而且如果你覺得頭部那里比較亂看不清的話,可以用鼠標拉一個框,它就放大了。
6.3 KMER統計表格
fastp對5個堿基長度的所有組合的出現次數進行了統計,然后把它放在了一張表格中,表格的每一個元素為深背景白字,背景越深,則表示重復次數越多。這樣,一眼望去,就可以發現有哪一些異常的信息。
image
[圖片上傳失敗...(image-4e31df-1576761117648)]
從上面的KMER表格中,我們可以發現,GGGGG的顏色特別深,從鼠標移上去之后顯示的信息中我們可以發現它的出現次數是平均次數的12.8倍,這是不正常的,因為GGGGG的正常倍數應該在1倍左右。幸好我們有fastp,它可以過濾掉這種polyG,讓數值較多地回歸正常。
6.4過表達序列(overrepresented sequence)
fastp的最新版本(v0.12)提供了overrepresetned
sequence的分析,而且不但提供了這些overrepresented sequence的序列個數和占比,還提供了他們在測序cycles中的分布情況,這有利于分析各種問題。具體示例如下圖所示:
7結語
好了,本次fastp的介紹就到此結束了,fastp軟件還在不斷更新中,目前每星期都有新功能開發出來,所以要想了解fastp軟件的最新動態,請關注該軟件的github項目地址https://github.com/OpenGene/fastp
總結
以上是生活随笔為你收集整理的linux质控命令,质控软件 fastp的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: ckfinder java 配置_JAV
- 下一篇: Linux中DNS服务器地址查询命令ns