本文轉(zhuǎn)載自知乎 nanopor技術(shù)專欄，修改了其中部分代碼。
單分子納米孔測(cè)序
https://poretools.readthedocs.io/en/latest/index.html#
三代測(cè)序技術(shù)
測(cè)序技術(shù)
2019年8月綜述 納米孔測(cè)序
液體活檢的納米孔測(cè)序：肺癌患者無細(xì)胞DNA拷貝數(shù)變異的分析
納米孔測(cè)序及其臨床應(yīng)用
一個(gè)人的全基因分析流程
使用的軟件
使用HDFView 查看Fast5文件格式。
https://www.hdfgroup.org/downloads/hdfview/
使用ont_fast5_api軟件對(duì)fast5文件進(jìn)行拆分與合并
ont_fast5_api 網(wǎng)址
https://github.com/nanoporetech/ont_fast5_api
安裝
直接使用pip安裝
pip install ont-fast5-api
自行編譯安裝
git clone https://github.com/nanoporetech/ont_fast5_api
cd ont_fast5_api
python3 setup.py install
使用案例
多條合并為一個(gè)文件
single_to_multi_fast5 -i fast5_files/ -s multi -n 4000 --recursive
一個(gè)文件拆分為多條
multi_to_single_fast5 -i multi/ -s single --recursive -t 6
1.堿基識(shí)別工具
DeepNano
poretools 官網(wǎng)
poretools安裝
2.序列比對(duì)工具
GraphMap
MarginAlign
3.從頭組裝工具
nanocorrect，首先利用graph- based，greedy partial order aligner方法進(jìn)行糾錯(cuò)，然后利用Celera Assembler將糾錯(cuò)后的reads進(jìn)行組裝，最后利用nanopolish對(duì)組裝結(jié)果進(jìn)行進(jìn)一步提升
4、單核苷酸變異檢測(cè)工具
PoreSeq 低通量下，仍然有高準(zhǔn)確率
Nanopolish
MarginAlign中的marginCaller模塊
5、共有序列的測(cè)序（consensus sequencing）方法
實(shí)戰(zhàn)
1.軟件安裝
比對(duì) minimap2 github
拼接 miniasm
gfatools
ngmlr
graphmap
拼接last
MarginAlign
過濾 filtlong
拼接 canu github
拼接 Smartdenovo
拼接 wtdbg2
拼接 NextDenovo
糾錯(cuò)、校正NextPolish
拼接 flye
基因組裝配 的質(zhì)量評(píng)估工具 quast
比對(duì) mummer 4
nanopack
pip install nanopack
基因組可視化 win+unix tablet
軟件給出官方地址，簡(jiǎn)單的自行安裝。比較復(fù)雜的會(huì)寫出詳細(xì)安裝過程
2. 質(zhì)控
mkdir nanoQC
#檢測(cè)測(cè)序質(zhì)量
nanoQC ../2.rawdata/minion/all..fastq.gz -o nanoQC
#統(tǒng)計(jì)質(zhì)量信息
NanoStat --fastq ../2.rawdata/minion/all.sra.fastq.gz --outdir statreports
2.1 質(zhì)控
NanoPlot --fastq ../2.rawdata/minion/all..fastq.gz -t 16 --maxlength 40000 --plots hex dot pauvre kde -o nanoplot
Nanoplot --summary sequencing_summary.txt --loglength -o summary
選項(xiàng)參數(shù)：
-t：線程數(shù)目
-o, --outdir：輸出結(jié)果目錄
-p, --prefix：輸出結(jié)果前綴
--color：點(diǎn)的顏色
--N50 表示在序列讀長(zhǎng)的直方圖中顯示N50的標(biāo)識(shí)
--title：標(biāo)題
--downsample ：在輸入文件中隨機(jī)抽取n條序列進(jìn)行處理
--minlength：忽略nbp以下的reads
-- fastq：輸入fastq格式文件
-f：圖片類型
--plots：繪圖類型，kde,hex,dot,pauvre
2.2 過濾（nanofilt和filtlong二選一即可）
使用nanofilt過濾，nanofilt無法識(shí)別壓縮文件，需要先解壓。
gunzip -c ../2.rawdata/minion/all.fastq.gz | NanoFilt -q 7 -l 1000 --headcrop 50 --tailcrop 50| gzip > clean.NanoFilt.fastq.gz
選項(xiàng)參數(shù)：
-l ,--length :過濾掉小于此長(zhǎng)度的序列
--maxlength ：過濾掉超過此長(zhǎng)度的序列
-q , --quality ：過濾掉低質(zhì)量序列
--minGC：過濾掉GC含量小于此百分比的序列
--maxGC：過濾掉GC含量大于此百分比的序列
--headcrop：從頭部切掉n bp
--tailcrop：尾部切掉n bp
結(jié)果解讀
輸入一個(gè)壓縮格式的fastq文件，結(jié)果軟件處理，已經(jīng)將長(zhǎng)度小于1000，同時(shí)整條序列平均質(zhì)量值小于Q7的過濾掉，同時(shí)將每條序列首位各剪掉50bp。這些剩余的序列則為“clean data”?？梢允褂肗anoPlot重新進(jìn)行一下質(zhì)控。
使用filtlong過濾
filtlong --min_length 1000 --min_mean_q 7 ../2.rawdata/minion/all.fastq.gz | gzip >clean.filtlong.fq.gz
filtlong還可以以二代的數(shù)據(jù)為參考進(jìn)行過濾。
選項(xiàng)參數(shù)：
--min_length ：最短長(zhǎng)度
--min_mean_q：平均Q值
--keep_percent：保留最好數(shù)據(jù)的百分比，80%，直接寫80，不能寫0.8
--target_bases：保留多少數(shù)據(jù)，單位為bp
2.3 比對(duì)
2.3.1 使用minimap2進(jìn)行比對(duì)
構(gòu)建基因組的索引文件
minimap2 -d ref.mmi ref.fa ＃索引
比對(duì)
minimap2 -a ref.mmi reads.fq > alignment.sam ＃對(duì)齊
#等價(jià)于
minimap2 -ax map-ont ref.mmi reads.fq >alignment.sam
常用選項(xiàng)參數(shù)
主要分成五大類，索引（Indexing），回帖（Mapping），比對(duì)(Alignment),輸入/輸出（Input/Output）,預(yù)設(shè)值（Preset）。
-x ：非常中要的一個(gè)選項(xiàng)，軟件預(yù)測(cè)的一些值，針對(duì)不同的數(shù)據(jù)選擇不同的值
map-pb/map-ont: pb或者ont數(shù)據(jù)與參考序列比對(duì)；
ava-pb/ava-ont: 尋找pd數(shù)據(jù)或者ont數(shù)據(jù)之間的overlap關(guān)系；
asm5/asm10/asm20: 拼接結(jié)果與參考序列進(jìn)行比對(duì)，適合~0.1/1/5% 序列分歧度；
splice: 長(zhǎng)reads的切割比對(duì)
sr: 短reads比對(duì)
-d :創(chuàng)建索引文件名
-a ：指定輸出格式為sa格式，默認(rèn)為PAF
-Q ：sam文件中不輸出堿基質(zhì)量
-R ：reads Group信息，與bwa比對(duì)中的-R一致
-t：線程數(shù)，默認(rèn)為3
2.3.2
2.4 sam轉(zhuǎn)bam,排序
#samtools處理
samtools sort -@ 8 -o bam -o s0137.sorted.bam s1037.sam
samtools index s0137.sorted.bam
samtools faidx ref.fq
2.5tablet可視化基因組
tablet支持Windows，linux.
下面是Windows上的操作
2.5.1 將準(zhǔn)備好的文件，至少四個(gè)，參考序列fasta格式及索引fai文件，比對(duì)并排序后的bam文件及索引bai文件；
ref.fq
ref.fai
ref.gff
s1037.sorted.bam
s1037.sorted.bam.bai
2.5.2、打開tablet，首先導(dǎo)入bam文件，然后導(dǎo)入fasta文件，索引文件無需導(dǎo)入，但必須與對(duì)應(yīng)文件在同一目錄下，（tablet對(duì)中文支持不好,路徑不要有中文）；
2.5.3 、可以通過鼠標(biāo)調(diào)整顯示模式；
2.5.4、tablet還支持導(dǎo)入gff格式結(jié)果，方便查看固定區(qū)域
2.6 重新拼接基因組（多個(gè)軟件分別拼接后，再?gòu)闹羞x擇最優(yōu)的）
2.6.1 canu
軟件安裝一定要安裝github的版本。
canu -d canu -p canu genomeSize=3g maxThreads=96 -nanopore-raw ../4.filter/clean.filtlong.fq.gz >canu.log
選項(xiàng)參數(shù)：
-p：輸出前綴
-d：輸出結(jié)果目錄
-nanopore-raw：輸入的為沒有糾錯(cuò)過得nanopore數(shù)據(jù)
-num_threads：CPU線程數(shù)
genomeSize：設(shè)置預(yù)估的基因組大小，這用于讓Canu估計(jì)測(cè)序深度；
maxThreads：設(shè)置運(yùn)行的最大線程數(shù)；
rawErrorRate：用來設(shè)置兩個(gè)未糾錯(cuò)read之間最大期望差異堿基數(shù)；
correctedErrorRate：則是設(shè)置糾錯(cuò)后read之間最大期望差異堿基數(shù)，這個(gè)參數(shù)需要在 組裝 時(shí)多次調(diào)整；
minReadLength：表示只使用大于閾值的序列
minOverlapLength：表示Overlap的最小長(zhǎng)度。提高minReadLength可以提高運(yùn)行速度，增加minOverlapLength可以降低假陽(yáng)性的overlap。
總結(jié)
1、有糾錯(cuò)步驟；
2、部分基因組拼接效果比較好；
3、默認(rèn)會(huì)占用所有CPU，非常耗時(shí)，慎重使用；
4、有些數(shù)據(jù)無法拼接出結(jié)果。
2.6.2 miniasm拼接，gfatools轉(zhuǎn)換格式
使用miniasm拼接首先需要使用minim2將測(cè)序數(shù)據(jù)進(jìn)行自身比對(duì)，查找共有區(qū)域，生成paf格式文件。注意使用minimap2比對(duì)的時(shí)候一定要正確設(shè)置好-x選項(xiàng)，nanopore拼接需要使用ava-ont選項(xiàng)。然后使用miniasm進(jìn)行拼接，miniasm拼接不會(huì)直接生成fasta序列，而是會(huì)生成gfa格式文件。最后使用gfatools提出拼接好的序列。
minimap2 -x ava-ont -t 12 ../4.filter/clean.filtlong.fq.gz ../4.filter/clean.filtlong.fq.gz | gzip -1 > reads.paf.gz
miniasm -f ../4.filter/clean.filtlong.fq.gz reads.paf.gz >reads.gfa
gfatools gfa2fa reads.gfa >miniasm.fa
注意：第一步，是用minimap2對(duì)測(cè)序的文件自身進(jìn)行比對(duì)。和前面的比對(duì)到基因組上不一樣。
總結(jié)
1、首先利用minimap2比對(duì)，運(yùn)行速度非?？?；
2、軟件不能直接生成fasta格式，需要使用gfatools生成；
3、沒有糾錯(cuò)過程，堿基錯(cuò)誤率較高，需要后續(xù)進(jìn)行糾錯(cuò)；
4、拼出來的基因組可能比真實(shí)基因組要小，有時(shí)候小很多，需要注意。
2.6.3 Smartdenovo拼接（阮玨 出品）
# 生成腳本
smartdenovo.pl -c 1 -t 8 -J 500 -k 16 -p smartdenovo ../4.filter/clean.filtlong.fq.gz >smartdenovo.mak
#運(yùn)行腳本
make -f smartdenovo.mak
選項(xiàng)參數(shù)：
-p STR 輸出結(jié)果前綴，默認(rèn)wtasm
-e STR overlap方法，zmo或者dmo
-t INT 線程數(shù)，默認(rèn)為8
-k INT overlap連接kmer大小
-J INT 最小read長(zhǎng)度
-c INT 是否生成一致性序列
2.6.4 wtdbg2拼接 （阮玨 出品）
快速運(yùn)行模式，直接使用wtdbg2.pl腳本
wtdbg2.pl -t 96 -x ont -g 3g -o wtdbg2 reads.fa.gz
選項(xiàng)參數(shù)：
-t 線程數(shù)量
-x 數(shù)據(jù)模式（nanopore的數(shù)據(jù)選擇ont）
-g 基因組大小（根據(jù)實(shí)際物種決定）
-o 輸出目錄
總結(jié)：
1、對(duì)于nanopore數(shù)據(jù)，wtdbg2可能會(huì)產(chǎn)生比真實(shí)基因組小的集合。
2、當(dāng)輸入fasta和fastq格式的多個(gè)文件時(shí)，請(qǐng)先輸入fastq，然后再輸入fasta。否則，程序?qū)o法在fastq中找到'>'，并在一次讀取后附加所有fastq。
2.6.5 NextDenovo拼接基因組
支持py2和py3
先用pip安裝依賴包Psutil和Drmaa (Only required by running under non-local system)
#生成一個(gè)配置文件，可以直接從軟件安裝目錄下拷貝
cp /ifs1/Software/biosoft/NextDenovo/doc/run.cfg ./
#生成一個(gè)reads列表，名為input.fofn
ls ../../filter/clean.filtlong.fq.gz > input.fofn
#將input.fofn添加到配置文件run.cfg的input_fofn關(guān)鍵字后面，默認(rèn)不用修改
input_fofn = ./input.fofn
#運(yùn)行軟件
nohup time nextDenovo run.cfg &
拼接結(jié)束之后，最終結(jié)果隱藏的比較深，在以下目錄中，可以配合nextpolish進(jìn)行糾錯(cuò)，得到更優(yōu)的結(jié)果。
01_rundir/03.ctg_graph/01.ctg_graph.sh.work/ctg_graph00/nextgraph.assembly.contig.fasta
2.6.6 flye拼接基因組
flye --nano-raw ../4.filter/clean.filtlong.fq.gz -g 3g -o output -t 48 -i 3
常用選項(xiàng)參數(shù)：（不習(xí)慣長(zhǎng)選項(xiàng)的可以直接使用短選項(xiàng)）
--pacbio-raw ：輸入原始pacbio數(shù)據(jù)
--pacbio-corr ：輸入糾錯(cuò)后的pacbio數(shù)據(jù)
--nano-raw：輸入原始nanopore數(shù)據(jù)
--nano-corr ：輸入糾錯(cuò)后的nanopore數(shù)據(jù)
--genome-size：預(yù)估基因組大小，用于評(píng)估覆蓋深度
--out-dir：輸出結(jié)果文件路徑
--threads：cpu線程數(shù)據(jù)
--iterations：糾錯(cuò)迭代次數(shù)
--min-overlap：最小overlap連接大小
--meta: 拼接宏基因組數(shù)據(jù)
--plasmids: 拼接質(zhì)粒數(shù)據(jù)
輸出結(jié)果
最后結(jié)果目錄中有三個(gè)文件比較重要。
1、assembly.fasta ：最終拼接得到的基因組序列，fasta格式。
2、assembly_graph.{gfa|gv} ：拼接過程中用到的repeat graph。
3、assembly_info.txt：拼接結(jié)果統(tǒng)計(jì)信息，也可以自己?jiǎn)为?dú)使用seqkit工具統(tǒng)計(jì)。
2.7用quast評(píng)估多個(gè)拼接結(jié)果
quast評(píng)估案例
軟件適合以下應(yīng)用場(chǎng)景：
1、得到不同軟件拼接的基因組序列，想比較一下哪個(gè)軟件拼接效果更好；
2、同一軟件，比較不同選項(xiàng)參數(shù)拼接的結(jié)果，哪個(gè)更好；
3、比較拼接得到的基因組，與參考序列的相似性。
quast.py -r ref.fa canu.fa miniasm.fa wtdbg2.cns.fa smartdenovo.fa -o quast
參數(shù)詳情
-o --output-dir 輸出結(jié)果目錄
-r 參考序列文件
-g --genes 參考序列基因坐標(biāo)，一般BED或者GFF格式文件
-m --min-contig 最小contig長(zhǎng)度，也就是小于這個(gè)閾值的不參與計(jì)算
-t --threads 使用線程數(shù)目，默認(rèn)使用四分之一的cpu
--help 列出全部選項(xiàng)參數(shù)，大部分情況下，默認(rèn)這些選項(xiàng)即可
2.8 組裝結(jié)果優(yōu)化（polish）
2.8.0 組裝優(yōu)化前后對(duì)比
使用Mummer軟件包中dnadiff工具將糾錯(cuò)前后的序列進(jìn)行基因組的比對(duì)。dnadiff會(huì)直接給出一個(gè)統(tǒng)計(jì)報(bào)告，里面會(huì)列出兩條序列之間差異的部分，然后我們可以使用mummerplot工具對(duì)統(tǒng)計(jì)進(jìn)行進(jìn)行粗略的可視化。
#拼接結(jié)果優(yōu)化前后比較
dnadiff ../before.fasta after.fasta
mummerplot --filter --png -p all out.delta
gnuplot all.gp
-i 輸入測(cè)序reads
-d 需要糾錯(cuò)的拼接結(jié)果
-o 輸出結(jié)果文件
-m 芯片類型
-t 并行計(jì)算
2.8.1 利用medaka組裝結(jié)果糾錯(cuò)
安裝方式 pip install medaka
medaka_consensus -i ../4.filter/clean.filtlong.fq.gz -d ../6.quast/miniasm.fa -o medaka_output -m r941_min_high -t 4 >medaka.log
-i 輸入測(cè)序reads
-d 需要糾錯(cuò)的拼接結(jié)果
-o 輸出結(jié)果文件
-m 芯片類型
-t 并行計(jì)算
最終生成的 consensus.fasta則為糾錯(cuò)后的結(jié)果。
calls_to_draft.bam
calls_to_draft.bam.bai
consensus.fasta
consensus.fasta.split/
consensus_probs.hdf
優(yōu)化前后比較
dnadiff ../6.quast/miniasm.fa medaka/consensus.fasta
mummerplot --filter --png -p all out.delta
gnuplot all.gp
2.8.2 pilon糾錯(cuò)
下載wget https://github.com/broadinstitute/pilon/releases/download/v1.23/pilon-1.23.jar，直接是二進(jìn)制文件。
#minimap2比對(duì)
minimap2 -d miniasm.fa.min miniasm.fa
minimap2 -ax map-ont miniasm.fa.min ../4.filter/clean.filtlong.fq.gz >miniasm.sam
#對(duì)bam進(jìn)行處理
samtools sort -@ 4 -O bam -o miniasm.sorted.bam miniasm.sam
samtools index miniasm.sorted.bam
#利用pilon進(jìn)行糾錯(cuò)
java -Xmx32G -jar /share/softwares/pilon/pilon-1.23.jar --genome miniasm.fa --fix all --changes --bam miniasm.sorted.bam --output all_pillon --outdir all_pillon --threads 6 --vcf 2> pilon.log
參數(shù)詳情
常用選項(xiàng)參數(shù)：
--genome提供輸入?yún)⒖蓟蚪M
--frags 表示輸入 < 1kb 的文庫(kù)BAM --jumps 輸入 > 1kb 的文庫(kù)BAM
-unpaired 輸入非配對(duì)的BAM。
--output表示輸出的前綴
--outdir表示輸出文件夾
--changes 列出發(fā)生變化的部分,以FASTA形式保存
--vcf 以VCF形式保存。
--fix 聲明對(duì)參考基因組做哪方面的改進(jìn)，包括“snps”,”indels”,”gaps”,”local”, 默認(rèn)是”all”
利用二代數(shù)據(jù)進(jìn)行糾錯(cuò)
#對(duì)拼接結(jié)果建立索引
bwa index draft.fasta
#illumina比對(duì)排序建索引
bwa mem -t 2 draft.fasta ../0.rawdata/SRX5341420_1.fastq.gz ../0.rawdata/SRX5341420_2.fastq.gz | samtools view - -Sb | samtools sort - -@ 14 -o illumina.sorted.bam
samtools index illumina.sorted.bam
#利用pilon進(jìn)行糾錯(cuò)
java -Xmx32G -jar /share/softwares/pilon/pilon-1.23.jar --genome draft.fasta --fix all --changes --frag --genome draft.fasta --fix all --changes --frags illumina.sorted.bam --output assembly_pillon --outdir assembly_pillon --threads 6 --vcf 2> pillon.log
結(jié)果解析
可能是因?yàn)閖ava寫的程序的緣故（不太確定），pilon的運(yùn)行速度比較慢。最終會(huì)生成*
*_pillon.changes：展示哪些位點(diǎn)經(jīng)過了糾錯(cuò)；
*_pillon.fasta ：最終糾錯(cuò)得到的結(jié)果；
*_pillon.vcf ：vcf格式展示糾錯(cuò)的位點(diǎn)。
