??現(xiàn)在三代測序越來越便宜了，很多小伙伴都對三代數(shù)據(jù)分析有著需求和好奇，這次直接一次性把三代流程給到大家，快上手操練起來吧。

1、軟件安裝

略（需要的話請留言或者加群）

2、質(zhì)控（2選1）

（1）nanoQC

mkdir?nanoQC#檢測測序質(zhì)量nanoQC??輸入數(shù)據(jù)路徑.fastq.gz?-o?nanoQC#統(tǒng)計質(zhì)量信息NanoStat?--fastq??輸入數(shù)據(jù)路徑.sra.fastq.gz?--outdir?statreports

（2）NanoPlot

#輸出質(zhì)量圖譜NanoPlot --fastq 輸入數(shù)據(jù)路徑.fastq.gz ?-t 16 --maxlength 40000 --plots hex dot pauvre kde -o nanoplot#統(tǒng)計質(zhì)量信息Nanoplot --summary sequencing_summary.txt --loglength -o summary

參數(shù)：

-t：線程數(shù)目

-o, --outdir：輸出結(jié)果目錄

-p, --prefix：輸出結(jié)果前綴

--color：點的顏色

--N50 表示在序列讀長的直方圖中顯示N50的標識

--title：標題

--downsample ：在輸入文件中隨機抽取n條序列進行處理

--minlength：忽略nbp以下的reads

-- fastq：輸入fastq格式文件

-f：圖片類型

--plots：繪圖類型，kde,hex,dot,pauvre

? ? ? ? ?

3、過濾（2選1）

這部分的參數(shù)選擇應該參考上一步的統(tǒng)計結(jié)果

（1）Nanofilt

#nanofilt無法識別壓縮文件，需要先解壓gunzip?-c?輸入數(shù)據(jù)路徑.fastq.gz?|?NanoFilt?-q?7?-l?1000?--headcrop?50?--tailcrop??50|?gzip?>?clean.NanoFilt.fastq.gz

參數(shù)：

-l ,--length :過濾掉小于此長度的序列

--maxlength ：過濾掉超過此長度的序列

-q , --quality ：過濾掉低質(zhì)量序列

--minGC：過濾掉GC含量小于此百分比的序列

--maxGC：過濾掉GC含量大于此百分比的序列

--headcrop：從頭部切掉n bp

--tailcrop：尾部切掉n bp

注：輸出結(jié)果需要用”>”導入文件中，否則只是輸出到屏幕上

? ? ? ? ?

（2）Filtlong

#fitlong無需解壓直接使用，此外還可以使用二代測序數(shù)據(jù)作為參考進行過濾filtlong --min_length 1000 --min_mean_q 7 輸入數(shù)據(jù)路徑.fastq.gz | ?gzip >clean.filtlong.fq.gz

參數(shù)：? ? ? ? ?
--min_length ：最短長度

--min_mean_q：平均Q值

--keep_percent：保留最好數(shù)據(jù)的百分比，80%，直接寫80，不能寫0.8

--target_bases：保留多少數(shù)據(jù)，單位為bp

注：輸出結(jié)果需要用”>”導入文件中，否則只是輸出到屏幕上

4、比對

（1）minimap2

#構(gòu)建索引minimap2?-d?ref.mmi?ref.fa#比對minimap2?-a?ref.mmi?reads.fq?>?alignment.sam#或者minimap2?-ax?map-ont?ref.mmi?reads.fq?>alignment.sam

優(yōu)點主要是：將nanopore測序數(shù)據(jù)比對到參考序列上就是整個nanopore數(shù)據(jù)分析的核心，有多種比對功能，除了支持nanopre數(shù)據(jù)之外，還支持pacbio數(shù)據(jù)。比對模式可以是reads與reads之前比對，reads與基因組比對，基因組與基因組比對，以及短序列與基因組的比對，不同的比對有不同的作用

參數(shù)：

-x ：非常中要的一個選項，軟件預測的一些值，針對不同的數(shù)據(jù)選擇不同的值

map-pb/map-ont: pb或者ont數(shù)據(jù)與參考序列比對；

ava-pb/ava-ont: 尋找pd數(shù)據(jù)或者ont數(shù)據(jù)之間的overlap關(guān)系；

asm5/asm10/asm20: 拼接結(jié)果與參考序列進行比對，適合~0.1/1/5% 序列分歧度；

splice: 長reads的切割比對

sr: 短reads比對

-d :創(chuàng)建索引文件名

-a ：指定輸出格式為sa格式，默認為PAF

-Q ：sam文件中不輸出堿基質(zhì)量

? ? ? ? ?

5、排序

#samtools處理samtools sort -@ 8 -o bam -o s0137.sorted.bam s1037.samsamtools index s0137.sorted.bamsamtools faidx ref.fq

? ? ? ? ?

6、拼接?

這個部分的軟件五花八門，以后考慮出個測評，此處簡單舉兩個栗子

（1）Canu

canu ?-p 輸出前綴 -d 輸出結(jié)果目錄 ?genomeSize=5g maxThreads=96 -nanopore-raw 原始的nanopore輸入數(shù)據(jù) > ?canu.log

注意：建議使用這個軟件時，數(shù)據(jù)先糾錯后，再組裝。

優(yōu)點主要是：cano的安裝特別容易，使用的最普遍，會自動進行糾錯、修剪、組裝。

（2）Unicycler

unicycler?-1?二代數(shù)據(jù)_1.fastq.gz?-2?二代數(shù)據(jù)_2.fastq.gz?-l?三代數(shù)據(jù).fastq.gz?-o?output_dir

優(yōu)點主要是：可以同時做二代和三代組裝外，對于三代組裝的操作也更為簡單，軟件內(nèi)部即可完成組裝、糾錯、成環(huán)等一系列步驟。現(xiàn)已獲得上千引用。

? ? ? ? ?

7.拼接質(zhì)量查看*

（1）Quast（此步驟可選）

quast.py?-r?ref.fa?canu.fa?miniasm.fa?wtdbg2.cns.fa?smartdenovo.fa?-o?quast

注：該工具常用于比較多個組裝軟件的拼接質(zhì)量

8.組裝結(jié)果優(yōu)化

（1）Medaka

medaka_consensus?-i?原始reads數(shù)據(jù).fastq.gz?-d?組裝后的assembly數(shù)據(jù).fasta?-o?medaka_result?-m?r941_min_high_g360?-v?medaka.vcf?-t?24?>?medaka.log

參數(shù)：

????????-i 輸入測序 reads

-d 需要糾錯的拼接結(jié)果

??? ????-o?輸出結(jié)果文件

-m 芯片類型,需要選好

-t 線程數(shù)

（2）Pilon（該軟件建議與minimap2聯(lián)合用）

使用二代數(shù)據(jù)對三代數(shù)據(jù)polish

#對拼接后結(jié)果建立索引mkdir pilonln -s 拼接后數(shù)據(jù)路徑/consensus.fasta assembly.fasta
#對bam進行處理samtools sort -@ 4 -O bam -o miniasm.sorted.bam miniasm.samsamtools index miniasm.sorted.bam
#利用pilon進行糾錯java -Xmx32G -jar 軟件路徑/pilon-1.23.jar --genome assembly.fasta --fix all --changes --bam miniasm.sorted.bam --output all_pillon --outdir all_pillon --threads 24 --vcf 2> pilon.log

參數(shù)：

--frags: 表示輸入的是1kb以內(nèi)的paired-end文庫，

--jumps表示 大于1k以上的mate pair文庫,

--bam輸入的bam文件

-vcf: 輸出一個vcf文件，包含每個堿基的信息

--fix: Pilon將會處理的內(nèi)容，基本上選snps和indels就夠了

--variant: 啟發(fā)式的變異檢測，等價于--vcf

--fix all,breaks, 如果是polish不要使用該選項

--minmq: 用于Pilon堆疊的read最低比對質(zhì)量，默認是0。

? ? ? ? ?

注意，pilon可以重復運行進行多次糾錯

只需將pilon結(jié)果再建立索引進行輸入就行

bwa index pilon.fasta

剩余代碼與之前一致，只需要更改文件名即可

? ? ? ? ?

參考：

https://www.jianshu.com/p/1053bf88f210

https://github.com/lh3/minimap2#general

https://lh3.github.io/minimap2/minimap2.html#10

https://doi.org/10.1101/169557

https://link.zhihu.com/?target=https%3A//github.com/marbl/canu

https://link.zhihu.com/?target=https%3A//canu.readthedocs.io/en/latest/tutorial.html

https://cloud.tencent.com/developer/article/2139549

? ? ? ? ?

隨著苗苗班的粉絲越來越多，有好幾個小伙伴都在問我有沒有培訓或者課程。想著我們的粉絲確實很多是剛剛到生信領(lǐng)域的小苗苗，這個需求我之前一直倒是沒注意到。我初步打算第一課就開簡單實用的全基因組分析實戰(zhàn)，大概就講周六周日兩個全天（講不完就再加時間）。主要內(nèi)容就是手把手帶著大家從測序原理，環(huán)境配置，到軟件安裝，數(shù)據(jù)實戰(zhàn)。不過我也不知道有多少人會報名，想征求一下大家的意見，有需求的就麻煩舉個手（留言或者加微信?。。?，如果能湊夠一個小班就開課。或者大家有什么別的想聽的主題，或者意見建議都歡迎提出哦~

?

(1)參加培訓、（2）商務(wù)合作、(3)付費分析、(4)進交流群

請加以下微信并備注目的。

這里是苗苗班，希望和大家共同交流進步，一起成長為參天有根樹??