【原】LongQC | 長讀長 Pacbio/Nanopore 測序數(shù)據(jù)質(zhì)控神器

生信藥丸 2021-11-03

展開全文

LongQC: A Quality Control Tool for Third Generation Sequencing Long Read Data

由于PacBio RS II ,、PacBio sequel 和 Nanopore測序質(zhì)量欠佳,，如果能對長序列進行質(zhì)控,，則也能對測序原始數(shù)據(jù)具有初步把握,。
推進LongQC的原因：1. 運行簡單 2. 結果展示美觀。

安裝

安裝依賴包

conda install h5pyconda install -c bioconda pysamconda install -c bioconda edlibconda install -c bioconda python-edlib

安裝longQC & minimap

cd /path/to/download/git clone https://github.com/yfukasawa/LongQC.gitcd LongQC/minimap2-coverage && makecd ..LongQCPath=`pwd`

由于部分python版本問題,，可能還需要額外安裝sklearn的python包,。

pip install sklearn

參數(shù)詳解

python $LongQCPath/longQC.py sampleqc -h

#usage: LongQC sampleqc [-h] -o OUT -x preset [-t] [-n NSAMPLE] [-s SUF] [-c TRIM] [--adapter_5 ADP5] [--adapter_3 ADP3] [-f] [-p NCPU] [-d] [-m MEM] [-i INDS] [-b]  input

#input 輸入待質(zhì)控數(shù)據(jù)的格式，如fasta、fastq或pbbam#-o OUT 輸出文件的路徑

#-x preset 輸入待質(zhì)控數(shù)據(jù)的測序類型,，以便自動提供數(shù)據(jù)的接頭和olvp參數(shù),，測序類型有:pb-rs2, pb-sequel, pb-hifi, ont-ligation, ont-rapid, ont-1dsq

#-t 輸入轉錄本、RNA或CDA序列文件（如果有）#-n NSAMPLE 樣本序列的數(shù)量,，默認為5000#-s SUF 輸出文件的前綴名字#-c TRIM 輸入去掉的接頭的保留路徑,，如果沒有輸入，接頭不會保留,。#--adapter_5 ADP5 輸入5'接頭#--adapter_3 ADP3 輸入3'接頭#-f 關閉一些sensitive的設置,，能運行更快但是準確性下降#-p NCPU 輸入線程數(shù)，默認4個線程,，線程數(shù)要求>=4#-d, --db 輸入minimap2比對的db庫

#-m MEM, --mem MEM  亞數(shù)據(jù)集的內(nèi)存限制,，指定為千兆字節(jié)（>0和<2)，默認是0.5千兆字節(jié)

#-i INDS, --index INDS 給出minimap2 (-I)的索引大小,，單位為bp,。在小內(nèi)存的服務器上運行時這個參數(shù)要減少。默認是4g,。

#-b, --short 對于非常短和高錯誤的reads(500bp),，開啟高靈敏度設置。

運行

采用hifi測序原始數(shù)據(jù)試一試

python $LongQCPath/longQC.py sampleqc -x pb-hifi -o out_dir input_reads.bam

結果展示

analysis/ #minimap比對結果figs/ #分析圖片logs/ #分析日志longqc_sdust.txt #QV的統(tǒng)計情況QC_vals_longQC_sampleqc.jsonweb_summary.html #可視化結果

結果解讀

質(zhì)控的結果概要

1. Sample name：測序樣品的名字,，可通過參數(shù)-s設置
2. Yield：堿基的總數(shù)量
3. Number of reads：reads數(shù)量
4. Q7 bases：堿基質(zhì)量值（QV）大于等于7的堿基的比例（由于是Pacbio 測序數(shù)據(jù),，所以并沒有提供QV）
5. Longest read：最長reads的長度
6. Estimated non-sense read fraction：評估出來的non-sense reads所占比例，應該少于30%,。比例高表明: 1. 測序數(shù)據(jù)有問題 2. 覆蓋度不高

接頭信息

用Edlib算法去查找數(shù)據(jù)中的接頭信息,。由于Pacbio數(shù)據(jù)已經(jīng)自動去掉接頭，所以這個部分對Pacbio數(shù)據(jù)默認是隱藏的,。

1. Number of trimmed reads：在reads末端中發(fā)現(xiàn)的可能是接頭（75%以上的一致性）的序列長度,；如果出乎意料地低或者沒有檢測到接頭，可能是在接頭連接步驟上有問題,。
2. Max seq identity：接頭序列與序列之間的最大的一致性值,。如果接頭仍然存在于數(shù)據(jù)集中，這個值應該相當高(90%),。
3. Average trimmed length：接頭序列在reads中的平均末端位置

Reads長度

（Nanopore 左,，PacBio 右）

這個圖展示是reads長度分布，三代測序數(shù)據(jù)會顯示單獨的高峰,。如果是轉錄組測序數(shù)據(jù),，會有嚴格的大小選擇的數(shù)據(jù)（如PacBio數(shù)據(jù)）會向右移動顯示更高的峰。

1. Mean read length：樣品的reads預測的長度

2. N50：N50值

Per-reads的質(zhì)量評估

（圖1 nanopore ,，圖2 PacBio）

上圖將顯示了per reads的QV分布,，y軸為平均QV值,，x軸是reads的長度集。
理想情況下,，短reads和長reads的QV分布應該相似,，中值(綠色橫線)應該高于7。
圖1是開發(fā)者用DASCRUBBERR軟件計算出QV后,，再進行質(zhì)控的結果,。
值得注意的是，圖2 由PacBio測序的序列是無法計算出QV,。主要是因為PacBio的測序文件質(zhì)量值是由"!"作占位符,，并未有實質(zhì)的QV意義。

Per read 覆蓋度

上圖為Per read 覆蓋度分布圖,。如果數(shù)據(jù)沒有問題,，應該是單峰（圖1）。如果是多峰值,，可以考慮是否基因組或者材料的雜合情況比較高（宏基因組數(shù)據(jù)除外）,。

LongQC使用GMM(用于基因組)或混合Gaussian 和lognorm 分布(用于轉錄組)自動檢測這樣的峰值,，以便從背景水平中區(qū)別真實的峰值,。

上圖為覆蓋度與reads長度的關系圖，用來檢查覆蓋范圍是否有意外波動,。在基因組測序數(shù)據(jù)中,，覆蓋范圍的波動應該在一定范圍內(nèi)(默認為3 σ)。如果超出這個范圍,，表明數(shù)據(jù)存在污染,，文庫質(zhì)量低，PacBio超載等問題,。
圖1表明存在微弱波動（可忽略）,，圖2是正常的分布圖。

non-sense reads

(Nanopore 左,，PacBio 右)

正常情況下,，Normal reads 比non-sense reads的顯示更高的值。
理想情況下,，正常reads的中值（橙色線）應該高于7,。

有3種情況：

1. 低coverage:當正常reads和non-sense reads的QV中值都在綠色范圍，說明數(shù)據(jù)低coverage,。因為當數(shù)據(jù)覆蓋度不夠高,，正常的reads會因為測序深度不夠而比對不上，導致被歸為non-sense reads,。

2. 噪音數(shù)據(jù)集：當正常reads和non-sense reads的中值都低于7（都在紅色區(qū)域）,，這就推斷出數(shù)據(jù)集的噪聲太大,，進一步的分析可能會受到嚴重影響。

3. 好的數(shù)據(jù)：當non-sense reads的中值低于7（紅色區(qū)域）,，正常reads的中值高于7（綠色區(qū)域）,，說明數(shù)據(jù)是準確的，如圖1,。而Pacbio測序數(shù)據(jù)是無法生成QV 圖,，且無法評估出Non-sense reads（圖2）。

注意：

這個部分顯示的中值覆蓋率可能比使用參考基因組進行質(zhì)控的結果要小,。由于將reads比對到未更正的容易出錯的整個測序數(shù)據(jù)集上時,，比對是不太敏感的，因此覆蓋率會受到這種不太敏感的結果的影響,。由于上述影響,，估計的基因組/轉錄組大小往往大于實際大小。一般來說,，更好的數(shù)據(jù)提供更好的尺寸估計,。

GC分布

上圖使用了兩種不同的方法計算GC含量，其中藍色是對end to end全長的reads進行計算,，而紅色是對分塊的亞數(shù)據(jù)集進行計算,。

藍色顯示更清晰的分布，因為使用更長的序列使結果受到更小的偏差,。由于測序平臺不同也會導致相同的數(shù)據(jù)的reads長度不同,，從而導致全長reads不同，而紅色的分塊亞數(shù)據(jù)集對長度已經(jīng)標準化,，所以對測序或樣本差異更不敏感,，更可靠。

Flanking region分析

可以用來檢查特定序列的連接,，如接頭,。如果沒有人工序列（如接頭），在兩個圖中峰值應該顯示在0的位置,。如果觀察到類似于接頭的序列,，則用紅色虛線標出平均長度。第一張圖表示有接頭,，第二張圖表示沒有接頭

序列完整性

這里通過DUST算法顯示序列的完整性分數(shù)

背景知識

三代測序的質(zhì)控問題

二代測序采用邊合成邊測序的方法,，以Illumina為例，進行橋式擴增后,，在Flow cell的固相表面上獲得上百萬條成簇分布的雙鏈待測片段,。然后在Flow cell中加入熒光標記的dNTP和酶，由引物起始開始合成子鏈,。由于dNTP存在 3'-OH 以疊氮基團RTG修飾,，這會阻礙子鏈延伸,，使得每個循環(huán)只能測得一個堿基。測完一個堿基后,， Flow cell 通入液體洗掉多余的dNTP和酶,，使用顯微鏡的激光掃描特征熒光信號，根據(jù)圖像的顏色即可識別該堿基,，再用化學試劑將疊氮基團與熒光基團洗去,。重復此過程即可識別全部堿基。在進行堿基識別的同時,，軟件還根據(jù)預先設定的規(guī)則對每一個堿基進行質(zhì)量評估,，為每一個堿基賦予一個質(zhì)量值。不同的二代測序平臺,，如果其測序原理不同,，則定義質(zhì)量值QV所取的參數(shù)也不同。

一個熒光信號被測序軟件判讀為某一種堿基,，存在一定的概率(Pe)發(fā)生判斷錯誤,。QV就是用來衡量Pe高低的一個定量指標，其定義為：
QV=-10 x lg Pe
根據(jù)這一定義可知,，
如果堿基識別的準確率是99%,，則錯誤率是1%，QV是20,；
如果堿基識別的準確率是99.9%,，則錯誤率是0.1%,，QV是30,；
如果堿基識別的準確率是99.99%，則錯誤率是0.01%,，QV是40
FASTQ格式的每第四行表示這條序列的質(zhì)量值,，用ACSII碼表示。所以二代測序的質(zhì)控軟件可以根據(jù)質(zhì)量值QV來評估測序數(shù)據(jù)的質(zhì)量,。

而三代測序由于其隨機分布的堿基錯誤率,，單堿基的準確性不能直接用來衡量數(shù)據(jù)的質(zhì)量，所以Pacbio Sequel測序平臺的下機數(shù)據(jù)不會提供單堿基質(zhì)量值QV,。此外,，三代測序雖然讀長長，但是錯誤率高（~15%）,，而且跟二代測序所用的技術原理不同,。二代測序的質(zhì)控軟件如FastQC（依賴QV評估測序數(shù)據(jù)）不適合三代測序的質(zhì)控。當三代測序數(shù)據(jù)不提供QV,，且沒有參考基因組時,，要怎么進行質(zhì)控呢,？

LongQC的介紹

LongQC質(zhì)控軟件可以解決以上的問題。LongQC不依賴QV和參考基因組進行質(zhì)控,，它依靠覆蓋度（Coverage)和non-sense reads對測序數(shù)據(jù)進行評估,。具體的原理和步驟如下：

1. LongQC首先將測序數(shù)據(jù)集隨機分成幾個亞數(shù)據(jù)集

2. 用minimap將亞數(shù)據(jù)集與整個測序數(shù)據(jù)集比對，計算overlap和coverage

3. 由于在測序文庫中核酸序列隨機分布,，當有足夠的coverage,reads跟整個測序數(shù)據(jù)集之間存在overlap.所以可以通過coverage來評估測序數(shù)據(jù),。

4. 當overlap少于2個時，就認為這個reads低coverage,，將這種reads定義為non-sense reads,，這些non-sense reads不能比對到整個測序數(shù)據(jù)集。

LongQC中的一些解釋

non-sense reads

non-sense reads是指沒有比對到整個測序數(shù)據(jù)集的reads,，這些reads在整個數(shù)據(jù)中的分布是任意且分散的,，開發(fā)者認為這些reads是測序技術或建庫中導致的錯誤reads，或者是污染的reads,。通過計算這些non-sense reads在整個測序數(shù)據(jù)集所占的比例,，可以評估測序數(shù)據(jù)的質(zhì)量。

堿基質(zhì)量值QV

由于Sequel測序數(shù)據(jù)沒有提供堿基質(zhì)量信息,，所以在longQC質(zhì)控結果中關于QV的評估都是0,。而Nanopore測序數(shù)據(jù)有Q-scores，Q-Scores是用來說明一個給定的堿基比它周圍的堿基更好還是更差的依據(jù),。開發(fā)者通過DASCRUBBER軟件對數(shù)據(jù)生成了QV,，再進行質(zhì)控后可以看到QV的解析結果。