對(duì)于二代測序,可以使用fastQC軟件對(duì)原始下機(jī)數(shù)據(jù)質(zhì)量進(jìn)行全面的統(tǒng)計(jì)及繪圖义起,多樣本時(shí)可進(jìn)一步結(jié)合使用multiQC滴须。但是由于二代測序(如 illumina 平臺(tái))獨(dú)特的特性馏艾,例如讀長長度一致诊杆,可能包含測序接頭(Adapter) 和PCR 引起的片段重復(fù)(Duplication)等蝗锥。因此妓灌,同是fastq格式文件轨蛤,fastQC并不適合用來處理三代測序數(shù)據(jù),包括Oxford Nanopore(ONT)測序平臺(tái)原始下機(jī)的fastq文件虫埂。
一祥山、ONT下機(jī)數(shù)據(jù)質(zhì)量
當(dāng)前ONT測序質(zhì)量雖然有很大的改善,但準(zhǔn)確性依然不及二代測序掉伏,例如illumina或者BGIseq等缝呕。2018-2019年主流芯片R9.4 準(zhǔn)確率對(duì)于2D reads為94%,1Dreads僅為86% ,如下圖Fig.2b所示(1)斧散。
現(xiàn)在最新的R10.4.1及配套試劑 (2)供常,使用Dorado v0.2.5進(jìn)行高準(zhǔn)確度堿基轉(zhuǎn)換(High accuracy ,HAC)原始數(shù)據(jù)平均堿基質(zhì)量可達(dá)Q20(99.0%)鸡捐,進(jìn)行超高準(zhǔn)確度堿基轉(zhuǎn)換(Super accuracy 栈暇,SUP)原始數(shù)據(jù)平均堿基質(zhì)量可達(dá)Q23(99.5%),如果進(jìn)行Duplex測序箍镜,原始數(shù)據(jù)平均堿基質(zhì)量可達(dá)Q30(99.9%)源祈。
當(dāng)前Nanopore測序下級(jí)數(shù)據(jù)仍然以Q7作為數(shù)據(jù)過濾的一個(gè)標(biāo)準(zhǔn)。PacBio三代測序平臺(tái)由于有CCS校準(zhǔn)模式(也有參數(shù)可以設(shè)置)色迂,一般來說HiFi數(shù)據(jù)個(gè)人覺得不用做質(zhì)控香缺,但對(duì)于現(xiàn)在的ONT數(shù)據(jù)個(gè)人建議做個(gè)質(zhì)控查看一下數(shù)據(jù)質(zhì)量,如果有需要可以做一下reads過濾歇僧。相信隨著ONT的測序準(zhǔn)確度越來越高以及Duplex測序方式图张,今后也可能弱化質(zhì)控過程。
二诈悍、演示ONT數(shù)據(jù) (fastq)
PRJNA809646:https://www.ebi.ac.uk/ena/browser/view/PRJNA809646
Whole-Genome Sequence analysis of Clinically Isolated Carbapenem Resistant Escherichia coli from Iran (IRANCOLI ONT, SRR21721846祸轮,SRR21721848)。
#下載 MinION sequencing: IRANCOLI ONT, SRR21721846
$ wget -c -t 0 ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR217/046/SRR21721846/SRR21721846_1.fastq.gz #數(shù)據(jù)大小521M
#下載 MinION sequencing: IRANCOLI ONT, SRR21721848
$ wget -c -t 0 ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR217/048/SRR21721848/SRR21721848_1.fastq.gz #數(shù)據(jù)大小543M
#fastq文件改名
mv SRR21721846_1.fastq.gz SRR21721846.fastq.gz
mv SRR21721848_1.fastq.gz SRR21721848.fastq.gz
三写隶、數(shù)據(jù)質(zhì)控軟件
直接上結(jié)論
NanoComp 用以對(duì)下機(jī)序列(reads)數(shù)據(jù)質(zhì)控倔撞,Chopper用以數(shù)據(jù)過濾,剪切和去污染序列足以慕趴。
1.NanoComp 包含了 NanoStat痪蝇,NanoPlot 和 NanoQC 的內(nèi)容鄙陡,而且多樣本之間的比較和整合性比較好,所以多樣本 可以直接使用 NanoComp 躏啰,單個(gè)樣本可以運(yùn)行Nanoplot 趁矾。
2.對(duì)下機(jī)原始序列進(jìn)行過濾,剪切和污染序列的去除可以使用Chopper 给僵。
3.如果運(yùn)行 NanoPlot 毫捣,則會(huì)產(chǎn)生 NanoStat.txt 文件,結(jié)果和單獨(dú)運(yùn)行 NanoStat 的結(jié)果一樣帝际。直接運(yùn)行Nanoplot就行 蔓同。
Wouter De Coster, Rosa Rademakers 實(shí)驗(yàn)室的博后蹲诀,編寫了一系列軟件對(duì)長度長測序原始下機(jī)數(shù)據(jù)進(jìn)行統(tǒng)計(jì)質(zhì)控和對(duì)污染序列進(jìn)行去除的軟件斑粱,也是現(xiàn)在常用的ONT質(zhì)控軟件。
NanoComp: comparing multiple runs 比較多組長度長測序和比對(duì)數(shù)據(jù)脯爪。
NanoStat: statistic summary report of reads or alignments 長度長測序和比對(duì)數(shù)據(jù)的統(tǒng)計(jì)報(bào)告则北。
NanoFilt: filtering and trimming of reads 過濾和修剪長度長序列。
NanoLyse: removing contaminant reads (e.g. lambda control DNA) from fastq 從fastq文件中去除污染序列(e.g. lambda對(duì)照/參照DNA序列)痕慢。
NanoPlot: Plotting tool for long read sequencing data and alignments. 對(duì)長度長測序和比對(duì)數(shù)據(jù)的統(tǒng)計(jì)報(bào)告進(jìn)行可視化尚揣。
以上所有軟件都存放在Nanopack的github里面,可以通過一行安裝命令進(jìn)行安裝:
#一鍵安裝ONT系列質(zhì)控軟件
$ pip install nanopack
# NanoComp NanoFilt NanoLyse NanoPlot NanoStat
Nanopack: https://github.com/wdecoster/nanopack
備注:
1. Nanostat 已經(jīng)被運(yùn)行速度更快的 Cramino 替代了:
CRAMINO: A tool for quick quality assessment of cram and bam files, intended for long read sequencing 對(duì)長度長序列的cram和bam文件進(jìn)行質(zhì)量評(píng)估掖举。
2. NanoFlit 和 NanoLyse 被運(yùn)行速度更快的 chopper 整合替代了:
Chopper: A rust implementation combining NanoLyse and NanoFilt into one faster tool for filtering, trimming, and removing contaminants快骗。
3.nanoQC需要單獨(dú)安裝。
下面我們介紹每一個(gè)軟件的具體安裝和使用方法:
1. NanoPlot
主頁網(wǎng)址::https://github.com/wdecoster/NanoPlot
- 軟件安裝
#pip install
$ pip install NanoPlot
$ pip install NanoPlot --upgrade
#conda install
$ conda install -c bioconda nanoplot
以前通過conda安裝過拇泛,在新買的云服務(wù)器上滨巴,兩種方法都報(bào)錯(cuò)。于是使用singularity的方法下載鏡像, 如果大家也有安裝問題可以嘗試使用singularity鏡像俺叭。
#從官方docker鏡像地址下載nanoplot的singularity鏡像
$ singularity pull nanoplot.sif docker://nanozoo/nanoplot
#使用
$ singularity exec nanoplot.sif NanoPlot
- 軟件使用
#官方例子
$ NanoPlot --summary sequencing_summary.txt --loglength -o summary-plots-log-transformed
$ NanoPlot -t 2 --fastq reads1.fastq.gz reads2.fastq.gz --maxlength 40000 --plots dot --legacy hex
$ NanoPlot -t 12 --color yellow --bam alignment1.bam alignment2.bam alignment3.bam --downsample 10000 -o bamplots_downsampled
#實(shí)際使用
$ NanoPlot -t 12 --N50 --dpi 300 \
--fastq samplename.fq.gz \
--title samplename \
-o samplename
$ singularity exec nanoplot.sif NanoPlot -t 12 --N50 --dpi 300 \
--fastq samplename.fq.gz \
--title samplename \
-o samplename \
# -N50 在片段長度N50處劃線
# --dpi 設(shè)置圖片dpi
# -t 線程數(shù)
#--title 所有表的標(biāo)題
- 演示數(shù)據(jù)
$ NanoPlot -t 24 --N50 --dpi 300 \
--fastq SRR21721846.fastq.gz \
--title SRR21721846 \
-o SRR21721846_NanoPlot
$ NanoPlot -t 24 --N50 --dpi 300 \
--fastq SRR21721848.fastq.gz \
--title SRR21721848 \
-o SRR21721848_NanoPlot
- 結(jié)果展示
$ cd SRR21721846_NanoPlot/
#結(jié)果文件:
NanoPlot-report.html 主要報(bào)告文件
LengthvsQualityScatterPlot_dot.html
LengthvsQualityScatterPlot_kde.html
NanoStats.txt
Yield_By_Length.html
Non_weightedHistogramReadlength.html
WeightedHistogramReadlength.html
Non_weightedLogTransformed_HistogramReadlength.html
Non_weightedLogTransformed_HistogramReadlength.png
WeightedLogTransformed_HistogramReadlength.html
NanoPlot_20231120_2259.log
每個(gè)結(jié)果目錄中都有NanoPlot-report.html文件恭取,用瀏覽器打開即可查看結(jié)果報(bào)告索引,報(bào)告中包括所有統(tǒng)計(jì)數(shù)據(jù)和交互圖熄守。
摘要統(tǒng)計(jì)Summary statistics
LengthvsQualityScatterPlot_dot
點(diǎn)狀展示每條read的長度和質(zhì)量的分布蜈垮,兩側(cè)加柱狀圖進(jìn)一步呈現(xiàn)長度和質(zhì)量的分布情況≡U眨可以看出長度分布在2K以下攒发,質(zhì)量集中在Q13.
LengthvsQualityScatterPlot_kde
核密度估計(jì)(Kernel Density Estimation,kde)展示每條read的長度和質(zhì)量的分布晋南,兩側(cè)加柱狀圖進(jìn)一步呈現(xiàn)長度和質(zhì)量的分布情況惠猿。可以看出長度分布在2K以下负间,質(zhì)量集中在Q13.
Nonweighted Histogram of read lengths 無加權(quán)長度分布直方圖
Nonweighted LogTransformed Histogram of read lengths 無加權(quán)及取對(duì)數(shù)后的長度分布直方圖
測序數(shù)量較大偶妖,且長度分布極不均勻且偏短姜凄,只在底部看到一條線,或一個(gè)峰趾访。此時(shí)就需要將數(shù)據(jù)進(jìn)行以10為底對(duì)數(shù)轉(zhuǎn)換再觀察态秧。
Weighted Histogram of read lengths 加權(quán)后的長度分布
Weighted LogTransformed Histogram of read lengths 加權(quán)及取對(duì)數(shù)后的長度分布
X軸為長度,Y軸是堿基數(shù)量扼鞋,更好地看出不同長度上的堿基數(shù)量分布申鱼。
加權(quán)(weighted)的含義: They're weighted by the number of nucleotides per bin. As such longer reads count 'heavier'. The height of the bar, therefore, does not represent the number of reads, but the number of nucleotides in that bin.
Yield by length 長度產(chǎn)出圖
X軸為長度,Y軸為產(chǎn)量的頻率云头。一般為越長越少捐友。
2. NanoQC
主頁網(wǎng)址:https://github.com/wdecoster/nanoQC
- 軟件安裝
#pip 安裝
$ pip install nanoQC
#conda安裝
$ conda install -c bioconda nanoqc
- 軟件使用
nanoQC [-h] [-v] [-o OUTDIR] fastq
positional arguments:
fastq Reads data in fastq.gz format.
optional arguments:
-h, --help show this help message and exit
-v, --version Print version and exit.
-o, --outdir OUTDIR Specify directory in which output has to be created.
-l, --minlen int Minimum length of reads to be included in the plots
This also controls the length plotted in the graphs
from the beginning and end of reads (length plotted = minlen / 2)
-l 參數(shù)制定最短的reads長度限制
- 演示數(shù)據(jù)
$ nanoQC SRR21721846.fastq.gz \
-o SRR21721846_NanoQC
$ nanoQC SRR21721848.fastq.gz \
-o SRR21721848_NanoQC
- 結(jié)果展示
cd SRR21721846_NanoQC/
#結(jié)果文件:
nanoQC.html
NanoQC.log
輸出結(jié)果包含log文件和html報(bào)告文件(主要看html):
報(bào)告中包含read長度、堿基含量和堿基質(zhì)量盘寡,個(gè)人覺得這個(gè)可以作為頭尾堿基過濾的一個(gè)參考楚殿。
3. NanoStat
主頁網(wǎng)址:https://github.com/wdecoster/nanostat
如果運(yùn)行NanoPlot,則會(huì)產(chǎn)生NanoStat.txt文件竿痰,結(jié)果和單獨(dú)運(yùn)行NanoStat的結(jié)果一樣。
- 軟件安裝
#pip 安裝
$ pip install nanostat
#conda安裝
$ conda install -c bioconda nanostat
- 軟件使用
#官方例子
$ NanoStat --fastq reads.fastq.gz --outdir statreports
$ NanoStat --bam alignment.bam alignment2.bam
- 演示數(shù)據(jù)
$ NanoStat --fastq SRR21721846.fastq.gz -o SRR21721846_NanoStat -n SRR21721846.txt
$ NanoStat --fastq SRR21721848.fastq.gz -o SRR21721848_NanoStat -n SRR21721848.txt
# -o 要和 -n 一起連用砌溺,才能生成 SRR21721848_NanoStat 文件夾以及文件夾里的 SRR21721848.txt 文檔(結(jié)果)
- 統(tǒng)計(jì)結(jié)果
$ cat SRR21721848.txt
#SRR21721848 統(tǒng)計(jì)結(jié)果
General summary:
Mean read length 平均read長度: 2,390.9
Mean read quality 平均堿基質(zhì)量 : 11.9
Median read length 長度中位數(shù) : 1,691.0
Median read quality 堿基質(zhì)量中位數(shù): 12.8
Number of reads 讀長總數(shù): 247,524.0
Read length N50 累計(jì)半總長的片段大小(N50) : 3,819.0
STDEV read length: 2,335.6
Total bases 總堿基數(shù): 591,806,555.0
Number, percentage and megabases of reads above quality cutoffs
>Q5: 247524 (100.0%) 591.8Mb
>Q7: 247524 (100.0%) 591.8Mb
>Q10: 229057 (92.5%) 553.8Mb
>Q12: 157978 (63.8%) 404.2Mb
>Q15: 29039 (11.7%) 87.7Mb
Top 5 highest mean basecall quality scores and their read lengths
1: 20.2 (199)
2: 19.6 (1538)
3: 19.5 (1930)
4: 19.4 (3240)
5: 19.0 (2429)
Top 5 longest reads and their mean basecall quality score
1: 37983 (15.3)
2: 35820 (12.8)
3: 33832 (11.2)
4: 33532 (12.4)
5: 33501 (10.2)
4. Cramino
主頁網(wǎng)址:https://github.com/wdecoster/cramino
- 軟件安裝
#conda安裝
$ conda install -c bioconda cramino
- 軟件使用
$ cramino --ubam samplename.bam -t 12
需要比對(duì)后或者未比對(duì)的bam文件影涉,這里就不做實(shí)際演示了,個(gè)人認(rèn)為跑一個(gè)NanoPlot就足以规伐。
- 結(jié)果展示
# Histogram for read lengths:
0-2000 ????????????????????????????????????????????????????????????????
2000-4000 ?????????????????????????????????
4000-6000 ????????????????????????????????
6000-8000 ???????????????????????????????
8000-10000 ???????????????????????????
10000-12000 ???????????????????????????
12000-14000 ??????????????????????????????
14000-16000 ???????????????????????????????????
16000-18000 ??????????????????????????????????????
18000-20000 ?????????????????????????????????????
20000-22000 ?????????????????????????????????
22000-24000 ???????????????????????????
24000-26000 ?????????????????????
26000-28000 ????????????????
28000-30000 ????????????
30000-32000 ?????????
32000-34000 ??????
34000-36000 ????
36000-38000 ??
38000-40000 ?
40000-42000 ?
42000-44000 ?
44000-46000
46000-48000
48000-50000
50000-52000
52000-54000
54000-56000
56000-58000
58000-60000
60000+
# Histogram for Phred-scaled accuracies:
Q0-1
Q1-2
Q2-3
Q3-4
Q4-5
Q5-6 ???
Q6-7 ????????????
Q7-8 ????????????????????????
Q8-9 ???????????????????????????
Q9-10 ???????????????????????
Q10-11 ?????????????????????????????
Q11-12 ??????????????????????????????????????????
Q12-13 ?????????????????????????????????????????????????????????
Q13-14 ??????????????????????????????????????????????????????????????????????
Q14-15 ???????????????????????????????????????????????????????????????????????????
Q15-16 ???????????????????????????????????????????????????????????????????
Q16-17 ???????????????????????????????????????????
Q17-18 ????????????????
Q18-19 ????
Q19-20 ?
Q20-21
Q21-22
Q22-23
Q23-24
Q24-25
Q25-26
Q26-27
Q27-28
Q28-29
Q29-30
Q30-31
Q31-32
Q32-33
Q33-34
Q34-35
Q35-36
Q36-37
Q37-38
Q38-39
Q39-40
Q40+
三蟹倾、數(shù)據(jù)質(zhì)控多樣本整合
NanoComp
主頁網(wǎng)址:https://github.com/wdecoster/nanocomp
Compare multiple runs of long read sequencing data and alignments. Creates violin plots or box plots of length, quality and percent identity and creates dynamic, overlaying read length histograms and a cumulative yield plot.
- 軟件安裝
#pip安裝
$ pip install NanoComp
#This script is written for Python3.
- 軟件使用
#官方示例
$ NanoComp --bam alignment1.bam alignment2.bam alignment3.bam --outdir compare-runs
$ NanoComp --fastq reads1.fastq.gz reads2.fastq.gz reads3.fastq.gz reads4.fastq.gz --names run1 run2 run3 run4
#實(shí)際樣本
$ nohup NanoComp -t 24 -f pdf \
--fastq SRR21721846.fastq.gz SRR21721848.fastq.gz \
--names SRR21721846 SRR21721848 \
-o NanoComp &
# -f 圖片以pdf的格式輸出
# -t 運(yùn)行線程數(shù)
- 結(jié)果展示
每個(gè)結(jié)果目錄中都有 NanoComp-report.html 文件,用瀏覽器打開即可查看結(jié)果報(bào)告索引猖闪,報(bào)告中包括所有統(tǒng)計(jì)數(shù)據(jù)和交互圖鲜棠。內(nèi)容其實(shí)和NanoPlot產(chǎn)生的結(jié)果一樣,只不過有多樣本之間的比較培慌。
摘要統(tǒng)計(jì)Summary statistics
總read數(shù)豁陆,總堿基數(shù),讀長N50數(shù)據(jù):
Read讀長吵护,對(duì)數(shù)轉(zhuǎn)換的Read讀長盒音,堿基質(zhì)量分布:
長度分布直方圖,標(biāo)準(zhǔn)化的長度分布直方圖馅而,加權(quán)后的長度分布:
取對(duì)數(shù)后的長度分布祥诽, 標(biāo)準(zhǔn)化及取對(duì)數(shù)后的長度分布,加權(quán)及取對(duì)數(shù)后的長度分布:
四瓮恭、數(shù)據(jù)過濾軟件
chopper
官方網(wǎng)站: https://github.com/wdecoster/chopper
Rust implementation of NanoFilt+NanoLyse, both originally written in Python. This tool, intended for long read sequencing such as PacBio or ONT, filters and trims a fastq file.
Filtering is done on average read quality and minimal or maximal read length, and applying a headcrop (start of read) and tailcrop (end of read) while printing the reads passing the filter.
基于平均Read質(zhì)量雄坪,最小最長Read讀長,頭尾堿基修剪
- 安裝方法
- 直接下載二進(jìn)制文件屯蹦,加入路徑维哈。
- Conda安裝盯漂。
#conda安裝
$ conda install -c bioconda chopper
- 使用方法
#使用說明
FLAGS:
-h, --help Prints help information
-V, --version Prints version information
OPTIONS:
--headcrop Trim N nucleotides from the start of a read [default: 0] 序列前端剪切
--maxlength Sets a maximum read length [default: 2147483647] 最大長度
-l, --minlength Sets a minimum read length [default: 1] 最小長度
-q, --quality Sets a minimum Phred average quality score [default: 0] 最小平均堿基質(zhì)量(Phred評(píng)分)
--tailcrop Trim N nucleotides from the end of a read [default: 0] 序列末端剪切
--threads Number of parallel threads to use [default: 4] cpu線程數(shù)
--contam Fasta file with reference to check potential contaminants against [default None] 去除污染序列
#使用示例
$ gunzip -c reads.fastq.gz | chopper -q 10 -l 500 | gzip > filtered_reads.fastq.gz
-
實(shí)際使用
根據(jù)質(zhì)控結(jié)果和實(shí)際需求調(diào)整
$ gunzip -c SRR21721846.fastq.gz | chopper -q 15 -l 500 | gzip > SRR21721846_filtered_reads.fastq.gz
參考文獻(xiàn):
- Wang, Y., Zhao, Y., Bollas, A., Wang, Y., & Au, K. F. (2021). Nanopore sequencing technology, bioinformatics and applications. Nature Biotechnology
- Raw read accuracy
- NanoPlot:三代納米孔測序數(shù)據(jù)質(zhì)量評(píng)估 - 劉永鑫
- 基因組組裝---Nanopore數(shù)據(jù)評(píng)估(nanoqc和NanoPlot套件工具)